# Robust Nonparametric Stochastic Frontier Analysis

Peng Zheng  Nahom Worku  Marlena Bannick  Joseph Dielemann  Marcia Weaver
Christopher Murray  Aleksandr Aravkin

###### Abstract

Benchmarking tools, including stochastic frontier analysis (SFA), data
envelopment analysis (DEA), and its stochastic extension (StoNED) are core
tools in economics used to estimate an efficiency envelope and production
inefficiencies from data. The problem appears in a wide range of fields – for
example, in global health the frontier can quantify efficiency of
interventions and funding of health initiatives.

Despite their wide use, classic benchmarking approaches have key limitations
that preclude even wider applicability. Here we propose a robust non-
parametric stochastic frontier meta-analysis (SFMA) approach that fills these
gaps. First, we use flexible basis splines and shape constraints to model the
frontier function, so specifying a functional form of the frontier as in
classic SFA is no longer necessary. Second, the user can specify relative
errors on input datapoints, enabling population-level analyses. Third, we
develop a likelihood-based trimming strategy to robustify the approach to
outliers, which otherwise break available benchmarking methods.

We provide a custom optimization algorithm for fast and reliable performance.
We implement the approach and algorithm in an open source Python package
‘sfma’. Synthetic and real examples show the new capabilities of the method,
and are used to compare SFMA to state of the art benchmarking packages that
implement DEA, SFA, and StoNED.

###### keywords:

Stochastic frontier , trimming , splines , nonlinear optimization

††journal: Journal of Econometrics

[IHME]organization=Institute for Health Metrics and
Evaluation,addressline=University of Washington, city=Seattle, postcode=98105,
state=WA, country=USA

[HMS]organization=Department of Health Metrics Sciences,addressline=University
of Washington, city=Seatle, postcode=98105, state=WA, country=USA

[Amath]organization=Department of Applied Mathematics,addressline=University
of Washington, city=Seatle, postcode=98105, state=WA, country=USA

[Biostat]organization=Department of Biostatistics,addressline=University of
Washington, city=Seatle, postcode=98105, state=WA, country=USA

Benchmarking tools, including stochastic frontier analysis (SFA) and data
envelopment analysis (DEA) are used to estimate an efficiency envelope and
production inefficiencies from data. The problem appears in a wide range of
fields – in the selected graphic above, the the frontier can be used to
quantify how efficiently countries use their GDP to extend life expectancy.

We propose a robust non-parametric stochastic frontier meta-analysis (SFMA)
approach. We use flexible basis splines and shape constraints to model the
frontier function, allow the user to incorporate reported errors for input
datapoints, and develop a likelihood-based trimming strategy to robustify the
approach to outliers.

We provide a custom optimization algorithm for fast and reliable performance,
and distribute the resulting method in the open source Python `sfma` package.
Synthetic and real examples show SFMA is competitive compared to state of the
art benchmarking packages in R, Python, and Stata.

We present a new semi-parametric benchmarking approach that generalizes
stochastic frontier analysis (SFA) to allow flexible spline modeling of the
frontier.

The meta-analytic focus allows the SFMA model to incorporate and estimate
different statistical errors, including reported sampling error, as well as
unknown non-sampling error and inefficiencies.

SFMA also features trimming functionality to remove outliers, that helps to
obtain useful fits in challenging real-world cases.

The SFMA package is available for general public use through the Python
package `sfma`, available at https://github.com/ihmeuw-msca/sfma.

##  1 Introduction

A goal in many econometric analysis is to estimate the production function –
the maximum amount of output obtainable from a set of inputs to a production
process – and the associated firm-specific production inefficiencies. The
framework is highly versatile, with a myriand benchmarking applications in the
energy industry [15, 25], agriculture [7, 8], education [35, 36], and
manufacturing [11, 30]. In the global health context, the output may be a
measure of healthcare access, the inputs may be investments in healthcare, and
inefficiencies may be country-specific gaps between the ideal production
function and observations [26]. Current approaches to the problem can be
classified along the lines of envelopment through linear programming,
exemplified by DEA, its stochastic extension StoNED, and likelihood-based
approaches such as SFA.

The approach in this paper, dubbed SFMA, addresses three key aspects:

  * 1.

Nonparametric modeling of the frontier using splines

  * 2.

Detailed modeling of statistical errors, including reported sampling error, as
well as unknown non-sampling error and inefficiencies

  * 3.

Automated outlier detection and removal.

To provide a robust software tool, we implement a customized optimization
algorithm that exploits problem structure and performs a single unified
analysis that determines outliers, estimates error statistics, and infers the
spline coefficients of the frontier and other covariate multipliers of
interest.

The SFMA approach can strongly influence our understanding and interpretation
of the results, particularly for complex real datasets like those presented
here: GDP vs life expectancy (LE), physician density vs universal health
coverage (UHC), and nurse density vs health coverage. For these datasets, SFMA
obtains helpful results, while a cross-section of tools available across
Stata, R, and Python fall short: StoNED fails to converge, SFA and DEA return
results that run counter to conventional wisdom in the field, as described
below.

  * 1.

For LE vs. GDP, DEA shows no change in LE for different levels of GDP, while
the SFA ‘frontier’ stays far below the upper bound of the data, and reports a
very steep (almost vertical) increase in life expectancy at low levels of GDP
per capita. The SFMA analysis suggests that LE increases gradually with GDP,
albeit more at lowest levels of GDP than higher. That result is made possible
due to the outlier-removal functionality available to SFMA.

  * 2.

For health coverage, the DEA analysis suggests that physicians have no effect
on UHC above 10 physicians per 10,000 population, a clearly falsifiable
conclusion. The SFA and SMFA analyses more plausibly show some benefit of more
physicians; SFA gain fails to reach the upper limits of the data.

In the sections below, we summarize existing solutions, their limitations,
discuss the contributions of SFMA, and present a roadmap for the paper.

###  1.1 Stochastic Frontier Analysis (SFA)

SFA [1] was a groundbreaking approach to reformulate deterministic estimation
of the frontier production function into a statistical model with stochastic
variation. From a statistical perspective, SFA is a mixed effects regression
that includes both symmetric random effects as well as asymmetric
inefficiencies for the units of analysis (i.e. ‘firms’). The frontier function
typically lies above most of the observed outputs for each firm, and is
estimated using a maximum likelihood formulation for the observed data
assuming both one-sided inefficiencies and random error. The frontier
estimates the best output that a firm can hope to achieve given input, while
fitted values of the inefficiencies for each firm estimate the improvable gaps
in outputs given inputs. Absent observation error, the frontier majorizes all
observations.

SFA analysis has been a topic of considerable interest, with a focus on
inefficiency modeling. The seminal paper [1] used half-normal and exponential
random variables for inefficiencies, developed an iterative approach to fit
the likelihood, and noted the sensitivity of the approach to outliers. Closed
form solutions for inefficiencies were provided by [23]. [19] used the gamma
distribution for inefficiencies, while [14] and [20] modeled inefficiencies
using non-parametric kernel density and dirichlet processes, respectively. [5]
allow the inefficiencies to be a function of firm-specific variables and time.
The Student’s T-Half Normal [38], the truncated normal [4], and the truncated
skewed Laplace [29] distributions have also been used to model the
inefficiencies. [18] use simulation to compute MLE estimates in cases where a
closed form solution of the likelihood are not available. the Student’s T-Half
Normal [38], the truncated normal [4], and the truncated skewed Laplace [29]
distributions have also been used to model the inefficiencies.

###  1.2 Data Envelopment Analysis and StoNED

Data envelopment analysis (DEA) uses linear programming to majorize input-
output observations from all the firms [9]. The approach is simple and robust,
resulting in thousands of cited applications. The downside is that there is no
way to account for error or uncertainty with DEA, setting up obvious
limitations in many applications and motivating stochastic extensions.

Stochastic Non-smooth Envelopment of Data (StoNED) is a semi-parametric
stochastic method that attemps to address the limitations of DEA to account
for noise and uncertainty. StoNED is implemented using two stages: (1)
applying Convex Nonparametric Least Squares (CNLS) to get a monotonically
increasing concave shape that best fits the data [24], and (2) estimating
conditional expected values using the residuals from the CNLS. In the second
stage, inefficiencies and variances are obtained using the Method of Moments
(MOM) and Pseudo-likelihood methods (QLE) [24] under specific distributional
assumptions.

###  1.3 Limitations

Current SFA approaches require the modeler to supply nonlinear transforms for
inputs and outputs so that the frontier can be framed as a linear model.
Finding these transforms can be a challenge, and in fact in our numerical
experiments with real data we found that even when a natural transform exists
(such as log-log), it still very useful to fit a nonparametric frontier to the
transformed dataset. Second, statistical models for inefficiencies cannot
handle outliers - the entire frontier can be disrupted by one or a few points
that deviate from the general pattern. This requires modelers to pre-process
the data, even more than for classic analyses such as linear regression.
Finally, existing numerical schemes struggle to converge for large datasets,
due to the numerical ill-conditioning of the SFA likelihood function, even in
the classic setting where inefficiencies are modeled as half-normal or
exponential. SFA requires working with the log of the complementary error
function, posing numerical challenges that have gone unaddressed until now.

The DEA approach is more robust, casting the frontier estimation problem as a
linear program, and obtaining a piecewise-linear frontier. Conversely, this
approach is unable to account for noise in the data, motivating stochastic
extensions, such as StoNED, which provides significant modeling flexiblity.
However, because StoNED uses a two-stage approach, the shape of the frontier
is fixed before distribution of the errors is accounted for, leaving only the
height of the frontier to be determined by residual analysis. In particular,
we cannot incorporate any statistics into the fitting process, since they are
unknown during the first stage. Despite the focus on stability in the
envelopment literature, we encountered difficulties using StoNED for the real
datasets of interest in the global health setting, as discussed in our
experiments.

Finally, neither SFA, nor DEA or StoNED are able to identify and remove
outliers using likelihood based methods.

###  1.4 Contributions

Here we develop a non-parametric extension to SFA that also optionally
incorporates reported standard errors, a common situation in meta-analysis. We
dub the method Stochastic Frontier Meta-Analysis (SFMA). We model the envelope
as a spline and estimate all coefficients and statistics in an all-at-once
likelihood-based approach. The interface allows the user to incorporate
reported standard errors to encode known sampling error, as well as allow for
unknown inefficiencies (positive random effects) and non-sampling error
(symmetric random effects). The method can be applied to general situations,
including population-level problems in global health.

The functionality of SFMA is compared to other methods in Table 1.

Table 1: Model comparison Model | Stochastic | Non-parametric | Specified Obs. SE | Outlier-Robust
---|---|---|---|---
DEA |  | ✓ |  |
SFA | ✓ |  |  |
StoNED | ✓ | ✓ |  |
SFMA | ✓ | ✓ | ✓ | ✓

##### Non-parametric frontier

We allow the frontier to be modeled using a spline. Scientists using the
stochastic frontier may have prior knowledge about the shape of the frontier
in relation to a set of inputs. While parametric models (e.g. linear or log-
linear) impose very strong assumptions, more typical prior knowledge concerns
the shape of the frontier (e.g. monotonically increasing or concave). We use
basis splines to capture nonlinear frontiers with optional constraints on
shape of the frontier. The main contrast to StoNED is that the shape is
estimated simultaneously with statistics of errors and inefficiencies using a
maximum likelihood approach.

##### Relative errors of input points

An analysis at the population level aggregates over input points with reported
specific errors, typically related to the sample sizes of input studies. This
information helps to distinguish more trustworthy points from larger
studies/locations from noisier points resulting from smaller sample sizes.
Meta-analysis uses these errors to test for the presence of non-sampling error
in the data, and to improve the accuracy of the aggregate. There is no way to
incorpprate reported uncertainty into SFA, DEA, or StoNED, but we can directly
use this information through the SFMA interface.

##### Robust Frontier Estimation Using Trimming

Real-world datasets often contain outlying firms. SFA, DEA, and StoNED do not
offer functionality for automatic detection of outliers. We develop a robust
extension of the least trimmed squares (LTS) method [33, 2, 39] for SFA that
makes it possible to estimate frontiers for datasets that may contain
outliers.

##### Customized Optimization Algorithm

We analyze the likelihood function required to fit SFA and SFMA, and highlight
a partial convex structure. We exploit this structure to develop a customized
block-coordinate descent algorithm, in effect separating a general nonconvex
ill-conditioned problem into a convex problem and a simple scalar nonconvex
problem. The specialized algorithm converges quickly on difficult problems,
whereas naive approaches using available tools (e.g. `scipy.minimize`) fail to
solve the SFA, or even evaluate the relevant likelihood function, on simple
examples.

###  1.5 Roadmap

The paper proceeds as follows. In Section 2 we discuss the statistical model
for SFMA, derive the likelihood, and discuss basis-splines with shape
constraints. We also develop the trimming concept to robustify SFMA. In
Section 3 we analyze the maximum likelihood objectives, present a result on
partial convexity, and develop an algorithm for solving the original problem
and the trimmed extension. Numerical results using simulated data are
presented in Section 4.1, including comparisons of SFMA to alternative
approaches on standard datasets and datasets containing outliers. Applications
to real datasets are presented in Section 4.2. All algorithms are implemented
and released in the `sfma` Python package111https://github.com/ihmeuw-
msca/sfma.

##  2 SFMA Formulation and Likelihoods

In this section, we develop the SFMA model that captures classic SFA, as well
as meta-analysis and meta-analytic extensions of SFA analysis.

###  2.1 Statistical Model

The SFMA model is specified as follows:

$\displaystyle y_{i}$ $\displaystyle=\left\langle
x_{i},\beta\right\rangle+u_{i}-v_{i}+\epsilon_{i}$ (1)
$\displaystyle\epsilon_{i}$ $\displaystyle\sim N(0,\sigma_{i}^{2}),$
$\displaystyle\sigma_{i}^{2}\mbox{ known }$ $\displaystyle u_{i}$
$\displaystyle\sim N(0,\gamma)$ $\displaystyle\quad\gamma\mbox{ unknown }$
$\displaystyle v_{i}$ $\displaystyle\sim HN(0,\eta)$
$\displaystyle\quad\eta\mbox{ unknown }$

where $y_{i}$ is the output of firm $i$, $\left\langle
x_{i},\beta\right\rangle$ encodes the frontier using firm-specific covariates,
$\epsilon_{i}$ are sampling errors with known variances $\sigma^{2}_{i}$,
$u_{i}$ are random effects, that is, non-sampling errors with unknown variance
$\gamma$, and $v_{i}$ are inefficiencies with unknown variance $\eta$. The
data-generating mechanism for model (1) is illustrated in Figure 1.

Table 2:  Model specifications using (1) Parameter Settings | Model type | References
---|---|---
{$\sigma_{i}^{2}=0$}, $\gamma$ free, $\eta$ free | Classic SFA | [1]
{$\sigma_{i}^{2}$} given, $\gamma$ free, $\eta=0$ | Meta-analysis | [10, 39]
{$\sigma_{i}^{2}$} given, $\gamma=0$, $\eta$ free | Meta-analytic frontier | This paper
{$\sigma_{i}^{2}$} given, $\gamma$ free, $\eta$ free | Meta-analytic frontier |
| with non-sampling error | This paper
Figure 1: Data generating process for stochastic frontier meta-analysis:
deviations from the frontier due to random effects, inefficiencies, and study
errors. Observed data points are shown in black at the bottom. Generated with
half-normal inefficiencies with $\eta=1$, $\gamma=0.25$, study-specific errors
drawn uniformly between 0 and 0.75, one covariate drawn uniformly between 0.5
to 1.5, and the frontier function $f(x)=\log(x)+10$.

Model (1) captures traditional forms of the SFA, as well as meta-analytic
extensions we introduce for population-level analysis, see Table 2. We proceed
now with developing the approach for (1), for simplicity restricting
inefficiencies $v_{i}$ to be half-Normal.

###  2.2 A Non-Parametric Frontier Using B-Splines

We use a flexible non-parametric spline representation rather than pre-
specified functional forms of the frontier. This approach allows us to
maintain flexibility and still incorporate a priori knowledge about the shape
of the frontier. Common assumptions about the frontier function include
monotonicity (more input results in larger output) and concavity (diminishing
returns with increasing inputs), and are encoded into the spline model using
linear constraints as discussed below.

Basis splines (B-splines) give the nonlinear frontier function a simple
representation as a linear combination of basis elements. Every spline
specification has at least two knots, defining the range of the domain. Using
degree $p$ polynomials and $k$ knots at locations by $t_{0},...,t_{k}$, we can
represent any frontier as the linear combination of $p+k$ spline basis
elements, with coefficients $\beta\in\mathbb{R}^{p+k}$ (see appendix section
7.2 for details):

$f_{\beta}(t)=\sum_{j=1}^{p+k}\beta_{j}s_{j}(t).$ (2)

Equation (2) represents arbitrary functions as linear combinations of the
basis elements, and is equivalent to equation (1) with $x_{i}$ rows from a
spline design matrix, see (21) in section 7.2.

We can specify the shape of the frontier using linear constraints to enforce
concavity, convexity, or monotonicity of the spline. These constraints are
imposed using derivatives: monotonicity is encoded using first derivatives,
while convexity/concavity is captured by second derivatives. For a basis
spline, derivatives of all orders are linear functions of spline coefficients,
and so all constraints can be imposed through the simple representation
$C\beta\leq c$ over a fixed grid, where $C\beta$ returns values of the
derivatives on the grid, and $c$ is the vector specifying required bounds for
these derivatives to achieve the shape constraints.

###  2.3 Likelihood

The likelihood corresponding to (1) is given by

$\displaystyle p(\beta,$
$\displaystyle\gamma,\eta,u,v|y)=p(\beta,\gamma,\eta|u,v,y)\cdot p(u,v|y)$ (3)
$\displaystyle\propto\prod_{i=1}^{n}\frac{1}{\sqrt{\sigma^{2}_{i}\gamma}}\exp\left(-\frac{(y_{i}-\left\langle
x_{i},\beta\right\rangle-
u_{i}+v_{i})^{2}}{2\sigma_{i}^{2}}-\frac{u_{i}^{2}}{2\gamma}\right)f(v_{i}|\eta)$
$\displaystyle:=\prod_{i=1}^{n}\mathcal{L}_{i}(\beta,\gamma,\eta,u,v).$

The model for $v_{i}$ is left unspecified in (3). The half-normal model

$f(v_{i}|\eta)=\begin{cases}\frac{\sqrt{2}}{\sqrt{\pi\eta}}\exp\left(-\frac{v_{i}^{2}}{2\eta}\right)&v_{i}\geq
0\\\ 0&v_{i}<0\end{cases}$ (4)

is widely used in practice [4],[28],[32],[21]. We stay with this choice
throughout the paper, and refer the reader to other models in the literature,
including Gamma [19], the Student’s T-Half Normal [38], the truncated normal
[4], and the truncated skewed Laplace [29] distributions.

When $v_{i}$ are half-normal (4), likelihood $p(\beta,\gamma,\eta,u,v|y)$ in
(3) is given by

$\displaystyle\prod_{i=1}^{n}\frac{1}{\sqrt{\sigma^{2}_{i}\gamma\eta}}\exp\left(-\frac{(y_{i}-\left\langle
x_{i},\beta\right\rangle-
u_{i}+v_{i})^{2}}{2\sigma_{i}^{2}}-\frac{u_{i}^{2}}{2\gamma}-\frac{v_{i}^{2}}{2\eta}\right)$
(5)

with support for each $v_{i}$ on the positive real line. Integrating out the
random effects, and taking the negative log of the resulting distribution, we
get a closed form solution of the marginal likelihood that only depends on
$(\beta,\eta,\gamma)$. To write this likelihood compactly, we first define
$\operatorname{erfc}$ to be the complementary error function

$\operatorname{erfc}(z)=1-\frac{2}{\sqrt{\pi}}\int_{0}^{z}\exp(-t^{2})dt.$ (6)

and also define

$r_{i}(\beta)=y_{i}-\left\langle x_{i},\beta\right\rangle,\quad
V_{i}(\gamma,\eta)=(\gamma+\eta+\sigma_{i}^{2}).$

Integrating $v_{i}$ and $u_{i}$ out of (5), and taking the negative logarithm,
we obtain the negative log likelihood used to fit $(\beta,\eta,\gamma)$:

$\displaystyle\ell_{hn}(\beta,\gamma,\eta|y)$
$\displaystyle=-\ln\left(\int_{\mathbb{R}^{m}_{+}}\int_{\mathbb{R}^{m}}p(\beta,\gamma,\eta,u,v|y)\,dudv\right)$
(7) $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\frac{r_{i}(\beta)^{2}}{2V_{i}(\gamma,\eta)}+\frac{1}{2}\ln
V_{i}-\ln\left(\operatorname{erfc}\left(\frac{\sqrt{\eta}r_{i}(\beta)}{\sqrt{2(\gamma+\sigma_{i}^{2})V_{i}(\gamma,\eta)}}\right)\right)$

The $\ln\operatorname{erfc}$ function in (7) presents a technical challenge to
algorithmic development, and is discussed in Section 3 along with other
details.

The SFMA model is fit by minimizing (7) subject to linear inequality
constraints that capture features such as monotonicity and convexity/concavity
(see Section 2.2) using the algorithm described in Section 3.2. The estimated
coefficients $\hat{\beta}$ give us the estimated frontier and any other
covariate multipliers of interest, while estimated variances
$(\hat{\gamma},\hat{\eta})$ are used to estimate individual non-sampling
errors $\hat{u}_{i}$ and inefficiences $\hat{v}_{i}$ as described in Section
3.4. In the next section, we discuss the trimming concept that modifies (7) to
make the proposed SFMA approach robust to outliers.

###  2.4 Extensions for Trimming

Outliers make it difficult to estimate a frontier using SFA [1]. For example,
a large outlier at low input can easily skew the entire frontier. Trimming
estimators from robust statistics [34] seek to simulatenously estimate a model
and automatically discover such outliers. Though trimming has a long history
(see e.g. [33]), it has recently been extended to general learning models
[AravkinDavis2020] and to meta-analysis [39]. Here we extend the trimming
concept to the SFMA likelihood in order to get an outlier-robust benchmarking
method.

Writing the negative log likelihood problem in (7) as a sum over individual
contributions over the dataset

$\min_{\theta}\sum_{i=1}^{n}\ell_{i}(\theta)\quad\mbox{s.t.}\quad C\theta\leq
c$ (8)

with $\theta:=(\beta,\gamma,\eta)$, we now considering auxiliary weights
$w_{i}$, one per datapoint, with constraints that each weight is in the
interval $[0,1]$ and all the weights sum up to the specified number of
inliers:

$\arg\min_{\theta,w}\sum_{i=1}^{n}w_{i}\ell_{i}(\theta),\quad\mbox{s.t.}\quad
C\theta\leq c,\quad 0\leq w_{i}\leq 1,\quad 1^{T}w=h.$ (9)

If $h=n$, then problem (8) requires that every $w_{i}=1$ and we reduce to the
original SFMA likelihood. For $h<n$, the weights are preferentially placed on
the $h$ “best fitting” observations (inliers), as judged by the likelihood at
a given $\theta$; the likelihood is in turn is affected by which points are
outliered in this manner. Since we optimize over both $\theta$ and $w$, the
joint problem is highly nonconvex [AravkinDavis2019], but both theoretical and
empirical results show the approach is stable and useful. Automated strategies
for choosing $h$ is an open problem, and in practice the analyst chooses a
common level across datasets, e.g. fitting 90% of the data to hedge the
quality of the results against potential outliers while still getting
information from the bulk of the data [27]. The method we use to estimate (8)
leverages the developments in [39] for trimmed meta-analysis, and is discussed
in Section 3.3.

##  3 Optimization Methods and Analysis

In this section we describe the optimization methods used to solve the SFMA
likelihood (5) as well as the corresponding trimmed extension (9). In order to
do this, we discuss technical challenges and theoretical properties of the
$\ln\operatorname{erfc}$ function, which plays an outsized role in (5). This
section proceeds as follows. In Section 3.1 we discuss
$\ln\operatorname{erfc}$ characterizes its convexity. This leads to a
customized block-coordinate descent method for (5) in Section 3.2, where the
problem is decomposed into an ill-conditioned convex problem with constraints
(solved using an interior point method) and two simple scalar nonconvex
problems that are easily solved. Finally, the algorithm for solving the
trimmed likelihood (9) discussed in Section 3.3.

###  3.1 The log of the complement to the Gaussian Error

The $\ln\operatorname{erfc}$ function

$\ln\operatorname{erfc}(z)=\log\left(1-\frac{2}{\sqrt{2\pi}}\int_{0}^{z}\exp(-t^{2})dt\right)$

plays a central role in the likelihood formulation and analysis for the SFMA
application, see (5). Simply evaluating the $\ln\operatorname{erfc}$ function
presents technical challenges, and current packages in Python do not evaluate
it correctly: there are no standard implementations of
$\ln\operatorname{erfc}$ in numpy/scipy in Python or in base R, and
specialized routines in Python and R return wrong values for large inputs. For
example, the function erfc in scipy.special.erfc and the erfc function in the
pracma package in R will both return exactly 0 for $x\geq 27$, and
$\ln\operatorname{erfc}(x)$ will correspondingly return $-\infty$ for these
values in both implementations. This is likely because inputs in this range
are not meaningful for simple statistical tests, but the issue plays a major
role in an optimization setting, where it causes optimization techniques to
fail for even moderately sized problems.

In order to fix the implementation of the basic function, we used the
asymptotic expansion for large $x$ and then created our own stable
implementation of $\ln\operatorname{erfc}$ that uses the asymptotic expansion
for large values. The asymptotic expansion of $\operatorname{erfc}(x)$ is
given by [22]:

$\operatorname{erfc}(x)=\frac{\exp(-x^{2})}{\sqrt{\pi}x}\left(1+\sum_{n=1}^{\infty}(-1)^{n}\frac{(2n-1)!!}{(2x^{2})^{n}}\right)$
(10)

We use the first two terms to approximate $\operatorname{erfc}(x)$ for large
$x$:

$\operatorname{erfc}(x)\approx\frac{\exp(-x^{2})}{\sqrt{\pi}x}\left(1-\frac{1}{2x^{2}}\right)$
(11)

And for $\ln\operatorname{erfc}(x)$, we have

$\ln(\operatorname{erfc}(x))\approx-x^{2}+\ln\left(1-\frac{1}{2x^{2}}\right)-\ln(\sqrt{\pi}x)$
(12)

For our Python implementation, we use standard routines when $x<25$ and (11)
for $x\geq 25$. This yields a stable implementation for the objective function
and makes it possible to optimize the required likelihoods.

Here, we establish concavity of $\ln\operatorname{erfc}$ and its approximation
(12). This allows us to show that the likelihood (5) is convex with respect to
$\beta$.

######  Theorem 1 (Concavity of $\ln\operatorname{erfc}$ and its
approximation).

The function $\ln\operatorname{erfc}$ is concave, and so is its approximation
(12) for $x\geq 25$.

The proof is given in the Appendix.

######  Corollary 1 (Partial convexity).

The marginal likelihood (7) is convex with respect to $\beta$.

###### Proof.

The likelihood (7) depend on $\beta$ through compositions of affine222Affine
functions are linear functions with an offset, e.g. $f(x)=Ax-b$. Compositions
of affine functions with convex functions are convex. functions $r_{i}(\beta)$
with convex functions: affine, quadratic, and $-\ln\operatorname{erfc}$. ∎

In the next section, we describe the block-coordinate descent approach that
takes advantage of the partial convexity structure in Corollary 1.

###  3.2 Block-coordinate Descent

The partial convexity of (5) with respect to $\beta$ makes it possible to
solve the problem using a block-coordinate descent algorithm, that alternates
full solves in blocks $(\beta,\gamma,\eta)$ until convergence. The largest
block is $\beta$, which may contain multiple parameters, including spline
coefficients and covariate multipliers, as well as constraints $C\beta\leq c$.
Convexity of the objective with respect to $\beta$ makes it possible to
optimize over $\beta$ efficiently using an interior point method described in
Section 7.4 of the Appendix. The remaining blocks are scalars, and we easily
solve these problem using root finding. The convergence of the block-
coordinate scheme here follows from classic results.

######  Theorem 2 (Convergence of Block-Coordinate Descent for SFA Marginal
Likelihood).

Suppose that the solution of (5) with respect to $\beta$ is unique for each
value of $\eta$ and $\gamma$, e.g. if the design matrix for the problem
$X^{T}X$ has full rank. Then the block coordinate scheme

$\displaystyle\beta^{k+1}$
$\displaystyle=\arg\min_{\beta}\ell(\beta,\gamma^{k},\eta^{k})$
$\displaystyle\gamma^{k+1}$
$\displaystyle=\arg\min_{\gamma}\ell(\beta^{k+1},\gamma,\eta)$
$\displaystyle\eta^{k+1}$
$\displaystyle=\arg\min_{\eta}\ell(\beta^{k+1},\gamma^{k+1},\eta)$

converges to a stationary point of (5).

###### Proof.

The convergence scheme satisfies conditions required by [37][Theorem 4.1 case
(c)]. Since there are only three blocks, as long as the convex problem in
$\beta$ has a unique solution, we can apply the theorem. ∎

In the next section, we show how to solve the extended likelihood (9),
leveraging the block-coordinate descent algorithm in this section.

###  3.3 Trimming Method

The algorithm for trimming uses variable projection [6, 17, 3]. Considering
the joint likelihood (9), and let $\theta=(\beta,\gamma,\eta)$. Define the
value function $v(w)$ and optimal value $\theta(w)$ by

$\displaystyle v(w)$
$\displaystyle=\min_{\theta}\sum_{i=1}^{m}w_{i}\ell_{i}(\theta)\quad\mbox{s.t.}\quad
C\theta\leq c$ (13) $\displaystyle\theta(w)$
$\displaystyle=\arg\min_{\theta}\sum_{i=1}^{m}w_{i}\ell_{i}(\theta)\quad\mbox{s.t.}\quad
C\theta\leq c.$

The term variable projection refers to partially minimizing over $\theta$. We
use the block-coordinate solver in Section 3.2 to solve the problem for each
fixed $w$, reducing the overall problem to value function optimization

$\min_{w}v(w)\quad\mbox{s.t.}\quad w\in\Delta_{h}:=\\{w:0\leq w_{i}\leq
1,\;1^{T}w=h\\}.$ (14)

where $v(w)$ is differentiable with derivative given by [39, Theorem 2]

$\nabla
v(w)=\nabla_{w}\left(\sum_{i=1}^{m}w_{i}\ell_{i}(\theta)\right)|_{\theta=\theta(w)}=\begin{bmatrix}\ell_{1}(\theta(w))\\\
\vdots\\\ \ell_{m}(\theta(w))\end{bmatrix}.$ (15)

The top level algorithm is then simply a projected gradient method where at
each iteration we compute the gradient and project the updated point onto the
set $\Delta_{h}$ defined in (14):

$w^{+}=\mbox{proj}_{\Delta_{h}}(w-\alpha\nabla v(w)).$

Each evaluation of $\nabla v$ requires a full minimization step over the
constrained weighted likelihood with respect to $\theta$ using our block-
coordinate descent method. The capped simplex $\Delta_{h}$ is a closed convex
set with a simple projection; proximal gradient with line search converges in
this case for any positive step-size $\alpha$ [AravkinDavis2019].

###  3.4 Estimating Firm-Specific Inefficiencies

Once fixed effects $\theta$ have been estimated, we want to obtain estimates
of both heterogeneity and inefficiency from the joint likelihood (5). Using
$r_{i}$ to denote the residuals

$r_{i}=y_{i}-\left\langle x_{i},\beta\right\rangle,$

we optimize

$\displaystyle\min_{u_{i},v_{i}\geq
0}\frac{(r_{i}-u_{i}+v_{i})^{2}}{2\sigma_{i}^{2}}+\frac{u_{i}^{2}}{2\gamma}+\frac{v_{i}^{2}}{2\eta}$
(16)

Partially minimizing with respect to $u_{i}$, we obtain

$\displaystyle u_{i}(v_{i})$
$\displaystyle=\frac{\gamma(r_{i}+v_{i})}{\sigma_{i}^{2}+\gamma}$ (17)

which shows that $\hat{u}_{i}$ is fully determined by $\hat{v}_{i}$. To solve
for the latter, we plug the estimate back into the expression of interest,
obtaining

$\min_{u_{i}}\frac{(r_{i}-u_{i}+v_{i})^{2}}{2\sigma_{i}^{2}}+\frac{u_{i}^{2}}{2\gamma}=\frac{1}{2(\sigma_{i}^{2}+\gamma)}(r_{i}+v_{i})^{2}.$
(18)

This reduces the problem to

$\displaystyle\min_{v_{i}\geq 0}$
$\displaystyle\frac{1}{2(\sigma_{i}^{2}+\gamma)}(r_{i}+v_{i})^{2}+\frac{v_{i}^{2}}{2\eta}$
(19) $\displaystyle\hat{v}_{i}$
$\displaystyle=\max\left(0,\frac{-\frac{1}{\sigma_{i}^{2}+\gamma}r_{i}}{\frac{1}{\sigma_{i}^{2}+\gamma}+\frac{1}{\eta}}\right)$
$\displaystyle=\max\left(0,-\frac{\eta
r_{i}}{\sigma_{i}^{2}+\eta+\gamma}\right)$

##  4 Numerical Experiments

In this section we present results for synthetic and real datasets comparing
SFMA to competing alternatives.

###  4.1 Synthetic Experiments

We created four different simulations and compared the results across SFMA,
DEA, StoNED, and SFA packages.

####  4.1.1 Simulation specifications

In all simulations, $x_{i}$ denote exposures of the independent variable,
$y_{i}$ denote the output of the firm $i$, $u_{i}$ is the firm-specific
inefficiency with variance $\sigma^{2}_{u}$ while $\epsilon_{i}$ is a random
error with standard error $\sigma^{2}_{\epsilon}$. For simulations 1-3, we
used

$\displaystyle x_{i}$ $\displaystyle=U(0,1)$ $\displaystyle y_{i}$
$\displaystyle=3+\log(x_{i}+0.2)+\epsilon_{i}-u_{i}$ $\displaystyle u_{i}$
$\displaystyle=HN(0,\sigma_{u}^{2}),\quad\epsilon_{i}=N(0,\sigma_{\epsilon}^{2})$

while for simulation 4, we used

$y_{i}=\begin{cases}3+\log(x_{i}+0.2)+\epsilon_{i}-u_{i}&\text{if inliers}\\\
10+\log(x_{i}+0.2)+\epsilon_{i}-u_{i}&\text{if outliers}\end{cases}$

with parameters as recorded in Table 3. The simulations become increasingly
difficult, with Sim 3 containing a mixture of reported standard errors, and
Sim 4 also containing outliers.

Table 3:  Simulation Specifications Sim | $\sigma_{u}^{2}$ | $\sigma_{\epsilon}^{2}$ | Outlier $\%$ | Data Points
---|---|---|---|---
1 | $1$ | $0.2$ | $0\%$ | 200
2 | $1$ | $\sqrt{0.2x_{i}}$ | $0\%$ | 200
3 | $0.5$ | $\begin{aligned} &0.05(67\%\text{ of data})\\\ &1.0(33\%\text{ of data)}\end{aligned}$ | 0% | 210
4 | $0.5$ | $\begin{aligned} &0.05(67\%\text{ of data})\\\ &1.0(33\%\text{ of data)}\end{aligned}$ | 12.5% | 210

####  4.1.2 Results for Simulated data

Table 4:  RMSE values for different models Sim | Sim 1 | Sim 2 | Sim 3 | Sim 4
---|---|---|---|---
DEA | 0.017290 | 0.33918 | 1.08809 | 62.54372
SFA | 0.05109 | 0.16333 | 0.01679 | 0.40832
StoNED (QLE) | 0.00360 | 0.01170 | 0.01551 | 5.81288
StoNED (MOM) | 0.00378 | 0.01179 | 0.01399 | 0.31752
SFMA | 0.00223 | 0.00399 | 0.00186 | 0.93045
R-SFMA | - | - | - | 0.00105

We compared multiple methods across the four simulations described above. For
each example, we ran a monte carlo simulations with 200 random realizations.
In each scenario and realization, we compared DEA, SFA, two versions of StoNED
(MOM and QLE), our method without trimming (SFMA), and our method with
trimming (R-SFMA). We also compared the frontiers of production functions
generated to the true generating function using root mean squared error
(RMSE), summarized in Table 4. The best-performing model in each sim is
highlighted in bold.

Figure 2: Results for Sim 1. SFMA and StoNED methods are close to the true
frontier, with SFA below and DEA above. RMSE for results across 200
realizations are in Table 4.

Results for Sim 1 (homoskedastic errors) are shown in Figure 2. The methods
are comparable and reasonably close to the true frontier, with the exception
of SFA, which deviates because of choice of functional form. SFMA obtains the
best results across multiple random realizations, as summarized in Table 4.

Figure 3: Results for Sim 2. SFMA is closest to true frontier, while StoNED
methods and SFA are below and DEA is significantly above. RMSE for results
across 200 realizations are in Table 4.

Results for Sim 2 (heteroscedastic errors) are shown in Figure 3. Point sizes
are inversely proportional to the reported SE of the data points. SFMA has an
advantage here, as it directly uses this information in the likelihood. This
is reflected in both the realization shown in the figure and in the overall
RMSE results in Table 4.

Figure 4: Results for Sim 3. SFMA significantly outperforms all other methods
as it is able to incorporate relative errors reported for all datapoints. RMSE
for results across 200 realizations are in Table 4.

Results for Sim 3 (mix of different SE levels) are shown in Figure 4. Point
sizes are inversely proportional to the reported SE of the data points. The
fact that different SE levels persist across the domain of interest
underscores how helpful it is to use this information. SFMA has a significant
advantage compared to all the other methods.

Figure 5: Results for Sim 4. SFMA with trimming (dubbed sfma$\\_$trim in the
plot) recovers the true frontier, while remaining methods deviate from truth
because of outliers added to the dataset. RMSE for results across 200
realizations are in Table 4.

Results obtained in the presence of outliers are shown in figure 5. Here, SFMA
with trimming offers a unique capability that is not available in any of the
other methods. We also see the potential for huge errors with straightforward
envelopment analysis from the large deviation of the DEA result.

Table 5:  List of Model Parameters Sim | No. knots | Degree | Trimming (for SFMA Trim) | Constraints
---|---|---|---|---
Sim 1 | 7 | 3 | 0 | Increasing, Concave
Sim 2 | 7 | 3 | 0 | Increasing, Concave
Sim 3 | 7 | 3 | 0 | Inreasing, Concave
Sim 4 | 7 | 3 | 12.5% | Increasing, Concave

Table 6 shows the list of parameters used to generate the plots.

###  4.2 Real datasets

In this section, we consider several challenging datasets used for
benchmarking in the global health setting. The first dataset focuses on life
expectancy (LE) as a function of GDP, while the second looks at a measure of
universal health coverage (UHC) as a function of physicians or nurses per
10,000 people.

####  4.2.1 World bank data

The World Bank dataset333https://data.worldbank.org/indicator/SP.DYN.LE00.IN
contains life expectancy from birth as a function of GDP across multiple
countries and years around the world. This dataset can be used to benchmark
life expectancy across 204 countries over the time period 1960-2021. A
frontier analysis alows us to evaluate the inefficiency for any country-year
at a given GDP per capita, helping to benchmark how effective countries are in
utilizing their GDP to improve life expectancy.

In most datasets, ‘outliers’ are a subjective notion. In this case, however,
multiple uploads from the World Bank gave us a unique opportunity to confirm a
set of outliers that was corrected by the World Bank. In the two uploads shown
in Figure 6, we see a set of country-years with high life expectancy and low
GDP (left panel) in the 2021 upload that were later corrected (right panel) in
the 2024 upload. We use the left dataset for the analysis to showcase the
robust trimming capability of SFMA.

Figure 6: Life expectancy as a function of GDP (on log scale). World bank
originally posted the data on the left panel, which includes a swath of
country-years far away from the bulk of the data. The current version of the
data, shown in the right panel, no longer contains the artifacts. We used the
left dataset to fit frontiers, to show the outlier-robust functionality of the
current method.

Figure 7 shows the frontiers built using SFMA with trimming, compared to DEA
and SFA. The spline feature of SFMA is very helpful here, since the frontier
is clearly nonlinear even after the classic log-log transform. The trimming
capability is key for this example – by allowing the model to identify and
exclude 3% of the data as outliers, SFMA gets a very informative fit, staying
ahead of the main data cloud while ignoring the spurious data in the top left
quadrant. SFA fits a curve essentially through the middle of the data cloud,
while DEA stays above both the main cloud and the outliers. The SFA result was
so surprising that we also fit the data using the Stata package `sfcross`, and
the result is also shown in Figure 7, dubbed ‘Stata SFA Unweighted’. We also
attempted to apply the StoNED approach on this example, but the CNLS fit in
the first stage of the method failed.

Figure 7: Life expectancy as a function of GDP, log scale (top) and linear
scale (bottom). DEA creates a flat line across the top, whie R SFA and Stata
SFA fit the middle of the data here. SFMA with trimming is able to get a much
more useful frontier by ignoring only 3% of the data.

####  4.2.2 UHC data

The UHC dataset was compiled by IHME444ghdx.healthdata.org/record/ihme-
data/gbd-2019-uhc-effective-coverage-index-1990-2019 in order to estimate
global progress towards universal health coverage (UHC) in 204 countries and
territories in 1990, 2010, and 2019. We use this dataset to fit a frontier of
the UHC index as a function of Physician density per 10,000 and Nurse density
per 10,000, also compiled by IHME555ghdx.healthdata.org/record/ihme-
data/gbd-2019-human-resources-health-1990-2019.

In fitting the frontiers, a few country-years with a low density of physicians
score unusually well with respect to UHC. These cases confuse standard
methods. The SFA approach fits through the middle of the data. In an effort to
obtain better frontier fits, we used the `sfcross` package in Stata. The
package allows weighting, and for UHC examples, we indeed found that inverse
variance weighting, using variances provided by the UHC dataset, produced a
more useful result in the `sfcross` package than no weighting, as is seen in
both linear and log scales in Figure 8 and Figure 9. Nonetheless, there is a
tendency for all SFA fits to go through the middle of the data, just as the
DEA result is guaranteed to go across the tops of the outlying datapoints,
flattening the estimate of the underlying relationship. SFMA with trimming is
able to get a useful fit to the majority of the data by automatically
detecting and excluding 3% of the data as outliers. We were not able to obtain
fits for the StoNED package on this dataset, as the CNLS routine failed to
converge. Further details relating to the real data examples are presented in
Section 5.

Figure 8: Universal healthcare access (UHC) as a function of physicians per
10,000 population, log scale (top) and linear scale (bottom). DEA fits above
all of the data. R SFA and default Stata SFA do poorly on this datset, fitting
through the middle of the data. Stata SFA with inverse variance weights gets a
more useful frontier. SFMA is able to obtain the closest fit to the main cloud
of data, also automatically trimming 5% of observations.

Figure 9: Universal healthcare access (UHC) as a function of nurses per 10,000 population, log-scale (top) and linear scale (bottom). DEA fits above all of the data. R SFA and both default and weighted Stata SFAs do poorly on this datset, fitting through the middle of the data. SFMA is able to obtain a closer fit to the cloud of data, also automatically trimming 5% of observations.  Table 6:  List of Model Parameters Dataset | No. knots | Degree | Trimming | Constraints
---|---|---|---|---
GDP - LE | 7 | 2 | 3% | Increasing, Concave
UHC-Physician | 7 | 2 | 5% | Increasing, Concave
UHC-Nurse | 7 | 2 | 5% | Increasing, Concave

##  5 Discussion

We have developed a new semi-parametric approach to frontier estimation,
dubbed SFMA. The approach generalizes previous approaches to SFA, allowing for
a semi-parametric frontier specification via splines with shape constraints, a
flexible likelihood that accounts for both estimates and reported uncertainty,
and can be used to automatically detect and remove outliers.

To facilitate easy access to the new method we developed and implemented a
customized optimization algorithm, which exploits the partially convex
structure of the marginal likelihood function.

Numerical results show the utility of the modeling innovations. The value
added of the innovations is clear as the simulations become more difficult,
and trimming in particular is indispensable in the real applications, where
ignoring small subsets of the data is essential to obtaining useful frontiers.

The applications to real datasets also demonstrate strengths of the modeling
innovations. SMFA estimates a positive relationship between LE and GDP, with a
steep increase from GDP per capita of $148 to $403 and more gradual from $2981
to $8103 than SFA. For example, in Vietnam, which is on the frontier from 2006
to 2011, LE increased from 74.3 in 2006 to 74.9 in 2011 when GDP per capita
increased from $784 to $1,525 . In Spain, which is on the frontier in 2016 and
2019, LE increased LE from 83.32 to 83.49 when GDP per capita increased from
$26,523. to $29,564. The functional form of SFA imposes a well-known, nearly
vertical relationship between LE and GDP at lowest levels of GDP, but accurate
data don’t exist in this range of GDP, and the relationship is well below the
frontier in the range where data exist. DEA estimates of the relationship as
nearly flat, which is also contrary to observed increases.

SMFA estimates a relatively steeper increase in UHC from 0 to 10 physicians
per 100,000 population, and more gradual increase above 20 physicians. Data on
the frontier in the lowest range of physicians per 100,000 population
represent countries such as Maldives with UHC score of 52.2 in 1990 and 0.9
physicians, 53.79 in 1992 and 1.05 physicians. Data on the frontier in the
middle range of physicians per 100,000 population represent countries such as
Iceland with a UHC score of 95.6 to 95.2 in 2011 and 2016, and 42.9 and 39
physicians, respectively. Three countries with 60 or more physicians per
100,000 population were well-below the frontier: Bermuda in 2009 and 2011,
Cuba from 2004 to 2017, and Georgia from 2006 to 2017.

The results suggest that there are other factors associated with increases in
UHC besides the number of physicians per capita. SFMA estimates a relatively
steeper increase in UHC from 0 to 10 nurses per 100,000 population, and more
gradual increase above 60 nurses. Data on the frontier in the lowest range of
nurses per 100,000 population represent countries such as Palestine with UHC
score of 60 in 1991 and 1992 and 5.91 and 6.04 nurses respectively, and 62 in
1998 with 7.35 nurses. Data on the frontier in the middle range of nurses per
100,000 population represent countries such as Tunisia with a UHC score of 74
in 2005 and 17.47 nurses, 77 in 2013 and 23.90 nurses, 78 in 2014 and 2015,
and 24.48, and 25.05 nurses, respectively. In the upper range, Luxemburg has a
UHC score of 94 in 2016 and 2017, and 88.3 and 87.7 nurses per 100,000
population, respectively.

The methods developed here can be extended to any statistical models for
inefficiencies that allows a closed form marginal likelihood, such as modeling
inefficiencies using the Laplace distribution. A closed form marginal
likelihood is leveraged to get the efficient algorithm that drives practical
performance, leaving extensions where the marginal likelihood is not available
in closed form to future work.

##  6 Acknowledgements

This work was funded by the Bill and Melinda Gates Foundation.

## References

  * [1] D. Aigner, C. K. Lovell, and P. Schmidt.  Formulation and estimation of stochastic frontier production function models.  Journal of econometrics, 6(1):21–37, 1977.
  * [2] A. Aravkin and D. Davis.  Trimmed statistical estimation via variance reduction.  Mathematics of Operations Research, 45(1):292–322, 2020.
  * [3] A. Y. Aravkin and T. Van Leeuwen.  Estimating nuisance parameters in inverse problems.  Inverse Problems, 28(11):115016, 2012.
  * [4] A. Baten, M. Rana, S. Das, and A. Khaleque.  Technical efficiency of some selected manufacturing industries in bangladesh: A stochastic frontier analysis.  Lahore Journal of Economics, 11(2):23 – 41, 2006.
  * [5] G. E. Battese and T. J. Coelli.  A model for technical inefficiency effects in a stochastic frontier production function for panel data.  Empirical economics, 20(2):325–332, 1995.
  * [6] B. M. Bell, J. V. Burke, and A. Schumitzky.  A relative weighting method for estimating parameters and variances in multiple data sets.  Computational statistics & data analysis, 22(2):119–135, 1996\.
  * [7] I. Benedetti, G. Branca, and R. Zucaro.  Evaluating input use efficiency in agriculture through a stochastic frontier production: An application on a case study in apulia (italy).  Journal of Cleaner Production, 236:117609, 2019.
  * [8] Z. Bibi, D. Khan, and I. u. Haq.  Technical and environmental efficiency of agriculture sector in south asia: A stochastic frontier analysis approach.  Environment, Development and Sustainability, 23:9260–9279, 2021\.
  * [9] P. Bogetoft and L. Otto.  Benchmarking with dea, sfa, and r, volume 157.  Springer Science & Business Media, 2010.
  * [10] M. Borenstein, L. V. Hedges, J. P. Higgins, and H. R. Rothstein.  Introduction to meta-analysis.  John Wiley & Sons, 2021.
  * [11] T. Charoenrat and C. Harvie.  The efficiency of smes in thai manufacturing: A stochastic frontier analysis.  Economic Modelling, 43:372–393, 2014.
  * [12] C. De Boor, C. De Boor, E.-U. Mathématicien, C. De Boor, and C. De Boor.  A practical guide to splines, volume 27.  springer-verlag New York, 1978.
  * [13] Y. Fan, Q. Li, and A. Weersink.  Semiparametric estimation of stochastic production frontier models.  Journal of Business & Economic Statistics, 14(4):460–468, 1996\.
  * [14] G. Feng, C. Wang, and X. Zhang.  Estimation of inefficiency in stochastic frontier models: A bayesian kernel approach.  Journal of Productivity Analysis, 51(1):1 – 19, 2019.
  * [15] M. Filippini and L. O. Sánchez.  Applications of the stochastic frontier approach in the analysis of energy issues.  Economics and Business Letters, 3(1):35–42, 2014.
  * [16] J. H. Friedman.  Multivariate adaptive regression splines.  The annals of statistics, pages 1–67, 1991.
  * [17] G. Golub and V. Pereyra.  Separable nonlinear least squares: the variable projection method and its applications.  Inverse problems, 19(2):R1, 2003.
  * [18] W. Greene.  Fixed and random effects in stochastic frontier models.  Journal of productivity analysis, 23(1):7–32, 2005.
  * [19] W. H. Greene.  A gamma-distributed stochastic frontier model.  Journal of econometrics, 46(1-2):141–163, 1990.
  * [20] J. Griffin and M. Steel.  Semiparametric bayesian inference for stochastic frontier models.  Journal of Econometrics, 123(1):121 – 152, 2004.
  * [21] A. Guerrini, G. Romano, and C. Leardini.  Economies of scale and density in the italian water industry: A stochastic frontier approach.  Utilities Policy, 52:103 – 111, 2018.
  * [22] J. K. Hunter.  Asymptotic analysis and singular perturbation theory.  Department of Mathematics, University of California at Davis, pages 1–3, 2004.
  * [23] J. Jondrow, C. K. Lovell, I. S. Materov, and P. Schmidt.  On the estimation of technical inefficiency in the stochastic frontier production function model.  Journal of econometrics, 19(2-3):233–238, 1982.
  * [24] T. Kuosmanen and M. Kortelainen.  Stochastic non-smooth envelopment of data: semi-parametric frontier estimation subject to shape constraints.  Journal of productivity analysis, 38:11–28, 2012.
  * [25] B. Lin and H. Long.  A stochastic frontier analysis of energy efficiency of china’s chemical industry.  Journal of Cleaner Production, 87:235–244, 2015.
  * [26] R. Lozano, N. Fullman, J. E. Mumford, M. Knight, C. M. Barthelemy, C. Abbafati, H. Abbastabar, F. Abd-Allah, M. Abdollahi, A. Abedi, et al.  Measuring universal health coverage based on an index of effective coverage of health services in 204 countries and territories, 1990–2019: a systematic analysis for the global burden of disease study 2019.  The Lancet, 396(10258):1250–1284, 2020.
  * [27] C. J. Murray, A. Y. Aravkin, P. Zheng, C. Abbafati, K. M. Abbas, M. Abbasi-Kangevari, F. Abd-Allah, A. Abdelalim, M. Abdollahi, I. Abdollahpour, et al.  Global burden of 87 risk factors in 204 countries and territories, 1990–2019: a systematic analysis for the global burden of disease study 2019\.  The Lancet, 396(10258):1223–1249, 2020.
  * [28] A. Musau, S. C. Kumbhakar, Ørjan Mydland, and G. Lien.  Determinants of allocative and technical inefficiency in stochastic frontier models: An analysis of norwegian electricity distribution firms.  European Journal of Operational Research, 288(3):983 – 991, 2021\.
  * [29] N. Nguyen, A. K. Gupta, and J. Wang.  Modeling technical efficiency using truncated skewed laplace distribution.  International Journal of Statistics and Economics, 15(3):80 – 99, 2014.
  * [30] A. Ö. Önder, E. R. DELIKTAS, and A. Lenger.  Efficiency in the manufacturing industry of selected provinces in turkey: a stochastic frontier analysis.  Emerging Markets Finance and Trade, 39(2):98–113, 2003.
  * [31] N. Pya and S. N. Wood.  Shape constrained additive models.  Statistics and Computing, 25(3):543–559, 2015.
  * [32] U. K. Raghunathan, C. L. Escalante, J. H. Dorfman, G. C. W. Ames, and J. E. Houston.  The effect of agriculture on repayment efficiency: A look at mfi borrowing groups.  Agricultural Economics, 42(4):465 – 474, 2011.
  * [33] P. J. Rousseeuw.  Least median of squares regression.  Journal of the American statistical association, 79(388):871–880, 1984.
  * [34] P. J. Rousseeuw and A. M. Leroy.  Robust regression and outlier detection, volume 589.  John wiley & sons, 2005.
  * [35] M. Sahnoun and C. Abdennadher.  Education and corruption: A stochastic frontier analysis: Evidence from developed and developing countries.  Journal of the Knowledge Economy, 11(3):968–981, 2020.
  * [36] M. A. Titus and K. Eagan.  Examining production efficiency in higher education: The utility of stochastic frontier analysis.  Higher education: Handbook of theory and research, pages 441–512, 2016.
  * [37] P. Tseng.  Convergence of a block coordinate descent method for nondifferentiable minimization.  Journal of optimization theory and applications, 109(3):475–494, 2001.
  * [38] P. Wheat, A. D. Stead, and W. H. Greene.  Robust stochastic frontier analysis: A student’s t-half normal model with application to highway maintenance costs in england.  Journal of Productivity Analysis, 51(1):21 – 38, 2019.
  * [39] P. Zheng, R. Barber, R. J. Sorensen, C. J. Murray, and A. Y. Aravkin.  Trimmed constrained mixed effects models: formulations and algorithms.  Journal of Computational and Graphical Statistics, 30(3):544–556, 2021.

##  7 Appendix

###  7.1 Derivative computations

In this section we compute gradient and hessian of the half-normal SFMA
likelihood. Recall the expression:

$\ell(\beta,\eta,\gamma)=\sum_{i=1}^{m}\frac{r_{i}^{2}}{2v_{i}}+\frac{1}{2}\ln(v_{i})-\ln\operatorname{erfc}(z_{i}r_{i})=\sum_{i=1}^{m}\ell_{i}(\beta,\eta,\gamma)$
(20)

where,

$r_{i}=y_{i}-\left\langle x_{i},\beta\right\rangle,\quad
v_{i}=\sigma_{i}^{2}+\eta+\gamma,\quad\tau_{i}=\sigma_{i}^{2}+\gamma,\quad
z_{i}=\sqrt{\frac{\eta}{2\tau_{i}v_{i}}}.$

We first compute all partial derivatives with respect to $\beta$, $\eta$ and
$\gamma$.

$\frac{\partial r_{i}}{\partial\beta}=-x_{i},\quad\frac{\partial
v_{i}}{\partial\eta}=1,\quad\frac{\partial
v_{i}}{\partial\gamma}=1,\quad\frac{\partial\tau_{i}}{\partial\gamma}=1.$

Next for $z_{i}$, which only depends on $\gamma$ and $\eta$, we have

$\frac{\partial
z_{i}}{\partial\eta}=\frac{z_{i}}{2}\left(\frac{1}{\eta}-\frac{1}{v_{i}}\right),\quad\frac{\partial
z_{i}}{\partial\gamma}=\frac{z_{i}}{2}\left(-\frac{1}{\tau_{i}}-\frac{1}{v_{i}}\right).$

And for the second order derivatives we have,

$\displaystyle\frac{\partial^{2}z_{i}}{\partial\eta^{2}}$
$\displaystyle=\frac{z_{i}}{4}\left(\frac{1}{\eta}-\frac{1}{v_{i}}\right)^{2}+\frac{z_{i}}{2}\left(-\frac{1}{\eta^{2}}+\frac{1}{v_{i}^{2}}\right),$
$\displaystyle\frac{\partial^{2}z_{i}}{\partial\eta\partial\gamma}$
$\displaystyle=\frac{z_{i}}{4}\left(\frac{1}{\eta}-\frac{1}{v_{i}}\right)\left(-\frac{1}{\tau_{i}}-\frac{1}{v_{i}}\right)+\frac{z_{i}}{2}\left(\frac{1}{v_{i}^{2}}\right),$
$\displaystyle\frac{\partial^{2}z_{i}}{\partial\gamma^{2}}$
$\displaystyle=\frac{z_{i}}{4}\left(-\frac{1}{\tau_{i}}-\frac{1}{v_{i}}\right)^{2}+\frac{z_{i}}{2}\left(\frac{1}{\tau_{i}^{2}}+\frac{1}{v_{i}^{2}}\right).$

The gradient and Hessian with respect to $\beta$ for a single term in the
negative log likelihood sum is given by

$\displaystyle\frac{\partial\ell_{i}(\beta,\eta,\gamma)}{\partial\beta}$
$\displaystyle=\left(-\frac{r_{i}}{v_{i}}+h^{\prime}(z_{i}r_{i})z_{i}\right)x_{i},$
$\displaystyle\frac{\partial^{2}\ell_{i}(\beta,\eta,\gamma)}{\partial\beta^{2}}$
$\displaystyle=\left(\frac{1}{v_{i}}-h^{\prime\prime}(z_{i}r_{i})z_{i}^{2}\right)x_{i}x_{i}^{\top}$

The first derivatives with repsect to $\eta$ and $\gamma$ are given by

$\displaystyle\frac{\partial\ell_{i}(\beta,\eta,\gamma)}{\partial\eta}$
$\displaystyle=-\frac{r_{i}^{2}}{2v_{i}^{2}}+\frac{1}{2v_{i}}-h^{\prime}(z_{i}r_{i})r_{i}\frac{\partial
z_{i}}{\partial\eta},$
$\displaystyle\frac{\partial\ell_{i}(\beta,\eta,\gamma)}{\partial\gamma}$
$\displaystyle=-\frac{r_{i}^{2}}{2v_{i}^{2}}+\frac{1}{2v_{i}}-h^{\prime}(z_{i}r_{i})r_{i}\frac{\partial
z_{i}}{\partial\gamma}.$

The second derivatives are given by

$\displaystyle\frac{\partial^{2}\ell_{i}(\beta,\eta,\gamma)}{\partial\eta^{2}}$
$\displaystyle=\frac{r_{i}^{2}}{v_{i}^{3}}-\frac{1}{2v_{i}^{2}}-h^{\prime\prime}(z_{i}r_{i})r_{i}^{2}\left(\frac{\partial
z_{i}}{\partial\eta}\right)^{2}-h^{\prime}(z_{i}r_{i})r_{i}\frac{\partial^{2}z_{i}}{\partial\eta^{2}},$
$\displaystyle\frac{\partial^{2}\ell_{i}(\beta,\eta,\gamma)}{\partial\eta\partial\gamma}$
$\displaystyle=\frac{r_{i}^{2}}{v_{i}^{3}}-\frac{1}{2v_{i}^{2}}-h^{\prime\prime}(z_{i}r_{i})r_{i}^{2}\frac{\partial
z_{i}}{\partial\eta}\frac{\partial
z_{i}}{\partial\gamma}-h^{\prime}(z_{i}r_{i})r_{i}\frac{\partial^{2}z_{i}}{\partial\eta\partial\gamma},$
$\displaystyle\frac{\partial^{2}\ell_{i}(\beta,\eta,\gamma)}{\partial\gamma^{2}}$
$\displaystyle=\frac{r_{i}^{2}}{v_{i}^{3}}-\frac{1}{2v_{i}^{2}}-h^{\prime\prime}(z_{i}r_{i})r_{i}^{2}\left(\frac{\partial
z_{i}}{\partial\gamma}\right)^{2}-h^{\prime}(z_{i}r_{i})r_{i}\frac{\partial^{2}z_{i}}{\partial\gamma^{2}}.$
$\displaystyle\frac{\partial^{2}\ell_{i}(\beta,\eta,\gamma)}{\partial\beta\partial\eta}$
$\displaystyle=\left(\frac{r_{i}}{v_{i}^{2}}+h^{\prime\prime}(z_{i}r_{i})z_{i}r_{i}\frac{\partial
z_{i}}{\partial\eta}+h^{\prime}(z_{i}r_{i})\frac{\partial
z_{i}}{\partial\eta}\right)x_{i},$
$\displaystyle\frac{\partial^{2}\ell_{i}(\beta,\eta,\gamma)}{\partial\beta\partial\gamma}$
$\displaystyle=\left(\frac{r_{i}}{v_{i}^{2}}+h^{\prime\prime}(z_{i}r_{i})z_{i}r_{i}\frac{\partial
z_{i}}{\partial\gamma}+h^{\prime}(z_{i}r_{i})\frac{\partial
z_{i}}{\partial\gamma}\right)x_{i}.$

###  7.2 Basis Splines

For general background on splines and spline regression see [12] and [16]. For
a spline with polynomial degree $p$ and $k$ knots placed at locations along
the domain $x$, $t_{0},...,t_{k}$, we need $p+k$ basis elements to construct
the spline basis, denoted $s^{p}_{j}$. These basis elements can be constructed
recursively.

The $j$th element over domain $c_{j}^{i}$ is denoted by $s_{j}^{i}$, with
$1\leq j\leq i+k$.

The spline basis is generated recursively: By default we assume
$c_{0}^{i}=\emptyset$, $c_{i+k+1}^{i}=\emptyset$ and $s_{0}^{i}(t)=0$,
$s_{i+k+1}^{i}(t)=0$, for all $i\geq 0$.

  * 1.

$i=0$,

$c_{j}^{0}=[t_{j-1},t_{j}),\quad s_{j}^{0}(t)=\delta_{c_{j}^{0}}(t),\quad
j=1,\ldots,k$

  * 2.

$i>0$,

$c_{j}^{i}=c_{j-1}^{i-1}\cup c_{j}^{i-1},\quad
s_{j}^{i}(t)=s_{j-1}^{i-1}(t)l(t;c_{j-1}^{i-1})+s_{j}^{i-1}(t)r(t;c_{j}^{i-1}),\quad
j=1,\ldots,i+k$

where

$\displaystyle l(t;c)$ $\displaystyle=\begin{cases}(t-\inf c)/(\sup c-\inf
c),&c\neq\emptyset\\\ 0,&c=\emptyset\end{cases},$ $\displaystyle r(t;c)$
$\displaystyle=\begin{cases}(t-\sup c)/(\inf c-\sup c),&c\neq\emptyset\\\
0,&c=\emptyset\end{cases}.$

as illustrated in Figure 10.

Figure 10: Generation of bspline basis elements (orders 0, 1, 2).

With knowledge of the degree and knots placement, we can construct a design
matrix $X$ based on an input vector $x$, where the $j^{th}$ column of the
design matrix is given by the expression

$\displaystyle X_{\cdot,j}=\begin{bmatrix}s_{j}^{p}(t_{0})\\\ \vdots\\\
s_{j}^{p}(t_{k})\end{bmatrix}.$ (21)

####  7.2.1 Linear Tails

We allow an option to make the first and/or last segment of the spline linear.
The prior capabilities can then be used to set a prior for the slope of this
segment to be $0$ (i.e. flat). The estimated spline is then a best fit to the
data, subject to this specification.

Figure 11: Spline extrapolation. Left: linear extrapolation. Right: nonlinear
extrapolation.

Given knots $t_{0},\ldots,t_{k}$, with $k\geq 3$, we want linear spline
functions in domains $[t_{0},t_{1})$ and $[t_{k-1},t_{k})$. We implement this
extension in two steps:

  * 1.

Generate elements $s_{j}^{p}(t)$ for the interval $[t_{1},t_{k-1})$ with knots
$t_{1},\ldots,t_{k-1}$ and degree $p\geq 1$.

  * 2.

Extend the basis linearly to $[t_{0},t_{1})$ and $[t_{k-1},t_{k})$.

The key observation is that only $s_{1},s_{2}$ have nonzero derivative at
$t_{1}$, and only $s_{k+p-3},s_{k+p-2}$ have nonzero derivative at $t_{k-1}$:

$\begin{matrix}{\displaystyle}s_{1}^{p}(t_{1})=1&\frac{d}{dt}s_{1}^{p}(t_{1})=-\frac{p}{t_{2}-t_{1}}&s_{k+p-2}(t_{k-1})=1&\frac{d}{dt}s_{k+p-2}(t_{k-1})=\frac{p}{t_{k-1}-t_{k-2}}\\\
s_{2}^{p}(t_{1})=0&\frac{d}{dt}s_{2}^{p}(t_{1})=\frac{p}{t_{2}-t_{1}}&s_{k+p-3}(t_{k-1})=0&\frac{d}{dt}s_{k+p-3}(t_{k-1})=-\frac{p}{t_{k-1}-t_{k-2}},\end{matrix}$

Therefore we extend our spline to the first and last piece by

$\displaystyle\bar{s}_{j}^{p}(t)$
$\displaystyle=\begin{cases}s_{j}^{p}(t_{1})+\frac{d}{dt}s_{j}^{p}(t_{1})(t-t_{1}),&t\in[t_{0},t_{1})\\\
s_{j}^{p}(t),&t\in[t_{1},t_{k-1})\\\
s_{j}^{p}(t_{k-1})+\frac{d}{dt}s_{j}^{p}(t_{k-1})(t-t_{k-1}),&t\in[t_{k-1},t_{k})\end{cases},$
$\displaystyle\frac{d}{dt}\bar{s}_{j}^{p}(t)$
$\displaystyle=\begin{cases}\frac{d}{dt}s_{j}^{p}(t_{1}),&t\in[t_{0},t_{1})\\\
\frac{d}{dt}s_{j}^{p}(t),&t\in[t_{1},t_{k-1})\\\
\frac{d}{dt}s_{j}^{p}(t_{k-1}),&t\in[t_{k-1},t_{k})\end{cases}.$

In general, using linear head and/or tail pieces to extrapolate outside the
original domain or interpolate in the data sparse region is far more stable
that using higher order polynomials, see Figure 11.

####  7.2.2 Shape Constraints

We can use constraints to enforce monotonicity, convexity, and concavity.
Monotonicity across the domain of interest follows from monotonicity of the
spline coefficients. This relationship is derived for particular basis
constructions by [12], and has been used in the literature to enforce shape
constraints [31]. Current approaches work around the natural inequality
constraints by using additional ‘exponentiated’ variables. Instead we impose
these constraints directly as described below.

Focusing just on $\alpha$, the relationship $\alpha_{1}\leq\alpha_{2}$ can be
written as $\alpha_{1}-\alpha_{2}\leq 0$. Stacking these inequality
constraints for each pair $(\alpha_{i},\alpha_{i+1})$ we can write all
constraints simultaneously as

$\underbrace{\begin{bmatrix}1&-1&0&\dots&0\\\ 0&1&-1&\dots&0\\\
\ddots&\ddots&\ddots&\ddots&\vdots\\\
0&\dots&\dots&1&-1\end{bmatrix}}_{\bm{C}}\begin{bmatrix}\alpha_{1}\\\
\alpha_{2}\\\ \alpha_{3}\\\ \vdots\\\
\alpha_{n}\end{bmatrix}\leq\begin{bmatrix}0\\\ 0\\\ \vdots\\\ 0\end{bmatrix}.$

These constraints are directly imposed through the optimization interface,
along with any lower- and upper-limit constraints on $\alpha$.

##### Convexity and Concavity

For any $\mathcal{C}^{2}$ (twice continuously differentiable) function
$f:\mathbb{R}\rightarrow\mathbb{R}$, convexity and concavity are captured by
the signs of the second derivative. Specifically, $f$ is convex if
$f^{\prime\prime}(t)\geq 0$ is everywhere, an concave if
$f^{\prime\prime}(t)\leq 0$ everywhere. We impose linear inequality
constraints on the expressions for $f^{\prime\prime}(t)$ over each interval.
We can therefore easily pick any of the eight shape combinations given in [31,
Table 1], as well as imposing any other constraints on $\alpha$ (including
bounds).

###  7.3 Proof of Theorem 1

We first establish the global concavity of $\ln\operatorname{erfc}$. For
convienience, we denote $h:=\ln\circ\operatorname{erfc}$. We have

$\displaystyle\operatorname{erfc}(x)$
$\displaystyle=1-\frac{2}{\sqrt{\pi}}\int_{0}^{x}\exp(-t^{2})\,\mathrm{d}t$
(22) $\displaystyle\operatorname{erfc}^{\prime}(x)$
$\displaystyle=-\frac{2}{\sqrt{\pi}}\exp(-x^{2})$
$\displaystyle\operatorname{erfc}^{\prime\prime}(x)$
$\displaystyle=\frac{4x}{\sqrt{\pi}}\exp(-x^{2})$

And we know,

$h(x)=\ln\operatorname{erfc}(x),\quad
h^{\prime}(x)=\frac{\operatorname{erfc}^{\prime}(x)}{\operatorname{erfc}(x)},\quad
h^{\prime\prime}(x)=\frac{\operatorname{erfc}^{\prime\prime}(x)\operatorname{erfc}(x)-\operatorname{erfc}^{\prime}(x)^{2}}{\operatorname{erfc}(x)^{2}}$

We first consider the case $x\leq 0$. From the formulas (22), we know that
$\operatorname{erfc}$ is a positive and decreasing, and
$\operatorname{erfc}^{\prime\prime}(x)\leq 0$ when $x\leq 0$. Therefore,

$h^{\prime\prime}(x)=\frac{\operatorname{erfc}^{\prime\prime}(x)\operatorname{erfc}(x)-\operatorname{erfc}^{\prime}(x)^{2}}{\operatorname{erfc}(x)^{2}}\leq
0,\quad\text{when }x\leq 0.$

Next, we consider the case when $x>0$. Here we use the asymptotic expansion of
$h$ at $\infty$,

$\operatorname{erfc}(x)=\frac{2}{\sqrt{\pi}}\exp(-x^{2})\left(\frac{1}{2x}+\sum_{k=1}^{\infty}\frac{(-1)^{k}(2k-1)!!}{2^{k+1}x^{2k+1}}\right)$

By grouping terms in the infinite sum in consecutive pairs, we see every pair
is non-positive, which means the entire infinite sum is less than or equal to
$0$. Therefore we have

$\operatorname{erfc}(x)\leq\frac{2}{\sqrt{\pi}}\frac{\exp(-x^{2})}{2x}\quad\Rightarrow\quad\operatorname{erfc}^{\prime\prime}(x)\operatorname{erfc}(x)-\operatorname{erfc}^{\prime}(x)^{2}\leq
0$

From this, we immediately have that $h^{\prime\prime}(x)\leq 0$, when $x>0$.
Combining with the previous case, we have $h^{\prime\prime}(x)\leq 0$ for all
$x\in\mathbb{R}$, and $h$ must be concave.

Next, the approximation for large values of $x$ is the sum of two terms:

$\left(-x^{2}-\ln(\sqrt{\pi}x)\right)+\ln\left(1-\frac{1}{2x^{2}}\right)$

The second derivative of the first term is given by $-2+\frac{1}{x^{2}}$ and
is negative for any $x>\frac{1}{\sqrt{2}}$. The second derivative of the
second term is

$\frac{-8x^{2}-4}{(2x^{2}-1)^{2}}+\frac{2}{x^{2}}\approx-\frac{2}{x^{2}}$

and is always negative. Thus the approximation is concave for
$x>\frac{1}{\sqrt{2}}$.

###  7.4 Interior Point Method for Smooth Convex Functions with Constraints

Consider the problem

$\min_{x}f(x)\quad\mathrm{s.t.}\quad\ Cx\leq c.$

where $f$ is twice continuously differentiable and convex. We introduce slack
variables $s$ such that

$Cx+s=c.$

The Lagrangian for the problem using the slack variables is given by

$\mathcal{L}(x,s,\lambda)=f(x)+\lambda^{\top}(Cx+s-c),$

and at the optimal point, complementarity slackness conditions yield
$\lambda_{i}s_{i}=0$ for all $i$. The Karush-Kuhn-Tucker (KKT) optimality
conditions are given by

$0=F\left(s,\lambda,x\right)=\begin{bmatrix}Cx+s-c\\\ \lambda\odot s\\\ \nabla
f(x)+C^{\top}\lambda\end{bmatrix}$

The interior point update is obtained by forcing the constraints to be
strictly satistfied, using a log-barrier method or, equivalently, modifying
the nonlinear complementarity portion of the KKT condition to be strict with
respect to a parameter $\mu$:

$F_{\mu}\left(s,\lambda,x\right)=\begin{bmatrix}Cx+s-c\\\ \lambda\odot
s-{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\mu 1}\\\
\nabla f(x)+C^{\top}\lambda\end{bmatrix}$

We then apply a damped Newton iteration to drive $F_{\mu}$ to $0$, with each
iteration given by

$\underbrace{\begin{bmatrix}I&0&C\\\ \Lambda&S&0\\\
0&C^{\top}&H(x)\end{bmatrix}}_{\nabla F_{\mu}}\begin{bmatrix}\Delta s\\\
\Delta\lambda\\\ \Delta x\end{bmatrix}=-\underbrace{\begin{bmatrix}Cx+s-c\\\
\lambda\odot s-\mu 1\\\ \nabla f(x)+C^{\top}\lambda\end{bmatrix}}_{F_{\mu}}$

where $H(x)$ is the hessian $\nabla^{2}f(x)$. The Newton update is
implementable whenever $H(x)$ is positive semidefinite (which is always true
for a twice differentiable convex $f$), and in fact we can obtain explicit
expressions of the Newton update for each coordinate:

$\displaystyle\Delta x$ $\displaystyle=(H(x)+C^{\top}\Lambda
S^{-1}C)^{-1}(-f_{3}+C^{\top}S^{-1}(f_{2}-\Lambda f_{1}))$
$\displaystyle\Delta s$ $\displaystyle=-f_{1}-C\Delta x$
$\displaystyle\Delta\lambda$ $\displaystyle=-S^{-1}(f_{2}+\Lambda\Delta s)$

For each iteration, we want to find the step size $\alpha$ that keeps the dual
and slack variables strictly positive, while also ensuring descent in the
optimality condition:

$\displaystyle\lambda^{+}(\alpha)>0$ $\displaystyle s^{+}(\alpha)>0$
$\displaystyle\|F_{\mu}^{+}(\alpha)\|_{\infty}\leq(1-\gamma\alpha)\|F_{\mu}\|_{\infty}$

Here we use $\gamma=0.01$, which makes the last inequality easy to satisfy for
relatively large steps. Once we find the step $\alpha$, we take

$\displaystyle s^{+}(\alpha)$ $\displaystyle=s+\alpha\Delta s$
$\displaystyle\lambda^{+}(\alpha)$ $\displaystyle=\lambda+\alpha\Delta\lambda$
$\displaystyle x^{+}(\alpha)$ $\displaystyle=x+\alpha\Delta x$ $\displaystyle
F_{\mu}^{+}$ $\displaystyle=F_{\mu}(s^{+},\lambda^{+},x^{+})$

To make the condition $F_{\mu}$ approach the optimality conditions $F$ of the
original problem, we drive the $\mu$ parameter to $0$ using a simple homotopy
technique:

$\mu^{+}=0.1\frac{\langle s^{+},\lambda^{+}\rangle}{\\#(s^{+})}.$

In other words, at each iteration, $\mu$ is taken to be a fraction of how far
complementarity conditions are away from $0$ on average. Once the (nonlinear)
complementary slackness conditions are satisfied, $F=0$ will hold, since the
affine equations in $F$ and $F_{\mu}$ are satisfied in every Newton update.

###  7.5 StoNED formulation

Consider the model

$\displaystyle y_{i}$ $\displaystyle=f(x_{i})-u_{i}-v_{i}$

where $u_{i}>0$ is the inefficiency term and $v_{i}$ is the noise term. We
assume $u_{i}$ has mean $\mu>0$ and variance $\sigma_{u}^{2}$, $v_{i}$ has
mean $0$ and variance $\sigma_{v}^{2}$. The first stage in StoNED estimates
the functional form of the frontier using CNLS:

$\displaystyle\min_{v,\alpha,\beta}\sum_{i=1}^{n}v_{i}^{2}$ $\displaystyle
y_{i}=\alpha_{i}+\beta_{i}x_{i}+v_{i}$
$\displaystyle\alpha_{i}+\beta_{i}x_{i}\leq\alpha_{h}+\beta_{h}x_{i},\forall
h,i=1,...,n$ $\displaystyle\beta_{i}\geq 0,\forall i=1,..,n$

Once $v_{i}^{2}$ is obtained, there are different choices of mthods to obtain
noise and inefficiency terms [24], as summarized below.

####  7.5.1 Method of Moments (MoM)

MoM assumes normality of sampling error and a half normal inefficiency term,
and uses skewness formulas to properties to disentangle the noise from the
inefficiency term [13]. Given the estimated second and third moment
$\hat{M_{2}}$ and $\hat{M_{3}}$of the CNLS residuals, the variance parameters
are given by [24].

$\displaystyle\hat{\sigma_{u}}=\sqrt[3]{\frac{\hat{M_{3}}}{\left[\frac{2}{\pi}\right]\left[1-\frac{4}{\pi}\right]}}$
$\displaystyle\hat{\sigma_{v}}=\sqrt{\hat{M_{2}}-[\frac{\pi-2}{\pi}}]\hat{\sigma_{u}}$

####  7.5.2 Pseudolikelihood Estimation (PSL):

Estimating the standard deviations $\sigma_{u}$ and $\sigma_{v}$ can also be
obtained using PSL, which was first coined by [13]. Here, parameters are set
as $\sigma=\sigma_{u}+\sigma_{v}$ , and obtained by maximizing log-likelihood
as a function of $\lambda$ , which is given as

$\displaystyle\log
L(\lambda)=-n\log\hat{\sigma}+\sum_{i=1}^{n}\log\phi[\frac{-\epsilon_{i}\lambda}{\hat{\sigma}}]-\frac{1}{2\hat{\sigma}^{2}}\sum_{i=1}^{n}\hat{\epsilon_{i}}^{2}$
$\displaystyle\hat{\varepsilon}_{i}=\hat{\upsilon}_{i}-\frac{{\sqrt{2}\lambda\hat{\sigma}}}{\pi(1+\lambda^{2})^{1/2}}$
$\displaystyle\hat{\sigma}=[\frac{\frac{1}{n}\sum_{j=1}^{n}\hat{\upsilon}^{2}}{[1-\frac{2\lambda^{2}}{\pi(1+\lambda)}]}]^{1/2}$

Here $\phi$ is the standard normal cumulative distribution function.

####  7.5.3 Estimation of inefficiency term

Once an estimator of $\sigma_{u}$ is available from MoM or PSL, the production
frontier is obtained by shifting the CNLS estimate by a factor of
$\hat{\sigma_{u}}\sqrt{\frac{2}{\pi}}$.
The conditional expected value of the inefficiency is given by

$\displaystyle\hat{E}(u_{i}\mid\hat{\varepsilon}_{i})=-\frac{\hat{\varepsilon}_{i}\hat{\sigma}_{u}^{2}}{\hat{\sigma}_{u}^{2}+\hat{\sigma}_{v}^{2}}+\frac{\hat{\sigma}_{u}^{2}\hat{\sigma}_{v}^{2}}{\hat{\sigma}_{u}^{2}+\hat{\sigma}_{v}^{2}}\left[\frac{\phi(\hat{\varepsilon}_{i}/\hat{\sigma}_{v}^{2})}{1-\Phi(\hat{\varepsilon}_{i}/\hat{\sigma}_{v}^{2})}\right]$

where
$\hat{\varepsilon}_{i}=\hat{\upsilon}_{i}-\hat{\sigma}_{u}\sqrt{\frac{2}{\pi}}$
is the estimator for the composite error term [24].


Magnetocaloric Effect in Two-Dimensional Diluted Ising Model: Appearance of
Frustrations in the Ground State  A.V. Shadrina,∗, V.A. Ulitkoa, Y.D. Panova
aUral Federal University Named after the First President of Russia B.N.
Yeltsin, Yekaterinburg, Russia

∗ e-mail<EMAIL_ADDRESS>
###### Abstract

The magnetocaloric effect in a two-dimensional Ising model is considered for
different ratios between parameters of inter-site repulsion of nonmagnetic
impurities and exchange coupling. Classical Monte Carlo simulations on a
square lattice show that in case of weak coupling and at sufficiently high
concentrations of nonmagnetic impurities the long-range ferromagnetic ordering
breaks down to give isolated spin clusters in the ground state of the system.
This leads to appearance of a paramagnetic response in the system at the zero
temperature and nonzero entropy of the ground state. The feasibility to detect
frustration of ground state using data on the magnetic entropy variation is
discussed.

##  1 Introduction

The magnetocaloric effect (MCE) consists in release or absorption of heat as a
result of variation of an external magnetic field applied to a material.
Initially, the effect was used to reach temperatures below 1 K, but, since
materials exhibiting MCE at near room temperature have been discovered,
magnetic cooling has become an area of very active research. The MCE in
frustrated and low-dimensional systems attracts a special interest [1, 2, 3,
4]. In Ising-type two-dimensional systems, the dependence of MCE on the
parameter of magnetoelastic interaction for a square lattice was investigated
in [5], while the effects of shape and size of geometrically frustrated Ising
spin clusters in a triangular lattice have on the MCE were addressed in [6].
In this work, we consider a two-dimensional Ising system with a fixed
concentration of mobile nonmagnetic charged impurities. The dilute Ising model
is one key model [7, 8] in the theory of magnetic systems with quenched or
annealed disorder as well as in the thermodynamic theory of binary alloys and
mixtures of classical and quantum liquids. To describe our system, we use the
pseudo-spin formalism ($S=1$) in which for a given lattice site the states
with pseudo-spin projections $S_{z}=\pm 1$ correspond to two magnetic states
with spin projections $s_{z}=\pm 1/2$, while the state with $S_{z}=0$
corresponds to a charged nonmagnetic impurity. Write the Hamiltonian as
follows:

$H=-\tilde{J}\sum\limits_{\langle ij\rangle}S_{zi}S_{zj}+V\sum\limits_{\langle
ij\rangle}P_{0i}P_{0j}-h\sum\limits_{i}S_{zi},$ (1)

where $S_{zi}$ is the $z$-projection of pseudo-spin operator on a site,
$P_{0i}=1-S_{zi}^{2}$ is the projection operator on state $S_{zi}=0$,
$\tilde{J}=Js^{2}$, $J$ is an exchange integral, $s=1/2$, $V$ is the inter-
site interaction between impurities, $h$ is an external magnetic field,
$\langle ij\rangle$ – are the nearest neighbors, and the sum is taken over all
sites of the twodimensional square lattice. The concentration n of charged
nonmagnetic impurities is fixed so that $nN=\sum\limits_{i}P_{0i}=const$ ,
where $N$ is the number of lattice sites.

The system described by Hamiltonian (1) exhibits two types of phase diagrams
that are exemplified in Fig. 1. An example with a strong coupling
($\tilde{J}>V$) is shown in the left panel (Fig. 1a). As the temperature is
lowered, two consecutive phase transitions occur in the system in the range of
$0<n<0.6$. The first one is a transition from a high-temperature disordered
state into ordered ferromagnetic (FM) one diluted with randomly distributed
charged impurities. At low temperatures, mobile charged impurities condense
into "drops". This means that in case of strong coupling the dilute FM phase
is unstable with respect to phase separation (PS) when the FM matrix expels
impurities in order to minimize surface free energy associated with them. Weak
coupling, i.e., when $\tilde{J}<V$, is exemplified in the right panel (Fig.
1b). With the impurity concentration increasing, the FM spin ordering
transforms into charge ordering (CO) of nonmagnetic impurities. Further, we
limit our discussion to the range of $0<n<0.5$, because for $n>0.5$ there is
no long-range order in a spin system with a fairly high critical temperature.

Fig. 1: Phase diagrams for the system with Hamiltonian (1): (a)
$V/\tilde{J}=0.4$ (strong coupling) and (b) $V/\tilde{J}=4$ (weak coupling).
The temperatures corresponding to maximum heat capacity values, as calculated
by a Monte Carlo method, are shown with filled circles. Dotted lines
schematically indicate boundaries of disordered (NO), ferromagnetic (FM), and
charge-ordered (CO) phases, as well as a phase separation region (PS).

The competition between magnetic and charge orders can be implemented in
systems subject to disproportionation or instability with respect to charge
transfer fluctuations [9], e.g., as in cuprates. Typically, with concentration
n of impurities increasing, the magnetic ordering is broken, the critical
temperature is lowered, and the temperature corresponding to maximum values
for MCE parameters, such as isothermal variation of magnetic entropy and
adiabatic temperature variation, is lower as well. Changes in the maximum MCE
temperature related to changes in the chemical composition have been well
demonstrated in [10, 11, 12]. However, the phase state of the system described
by Hamiltonian (1) also depends on the ratio between the constants of exchange
and densitydensity interactions, as was shown for a similar spinpseudo-spin
model [13, 14, 15, 16]. In this work, we compare the MCE parameters in cases
of strong and weak coupling; we also explore the possibility of detecting
frustration in ground state of the system with weak coupling using magnetic
entropy data.

##  2 Calculation procedure for the MCE parameters using Monte Carlo data

The adiabatic temperature change $\Delta T_{ad}$ and isothermal magnetic
entropy change $\Delta S_{M}$ are the key parameters characterizing MCE.
$\Delta S_{M}$ is calculated by the equation:

$\Delta
S_{M}(T,h)=\int\limits_{0}^{T}\dfrac{C(T^{\prime},h_{m})-C(T^{\prime},0)}{T^{\prime}}dT^{\prime},$
(2)

where $C$ is a heat capacity. Maxwell relationship $(\partial S/\partial
h)_{T}=(\partial M/\partial T)_{h}$ gives us another expression for $\Delta
S_{M}$ and an equation for $\Delta T_{ad}$:

$\Delta S_{M}(T,h_{m})=\int\limits_{0}^{h_{m}}\left(\frac{\partial
M(T,h)}{\partial T}\right)_{h}dh$ (3)

$\Delta
T_{ad}(T,h_{m})=\int\limits_{0}^{h_{m}}\frac{T}{C(T,h)}\left(\frac{\partial
M(T,h)}{\partial T}\right)_{h}dh,$ (4)

where $M$ is a magnetization and $h$ is an external magnetic field. In case of
second-order phase transitions, expressions (2) and (3) must lead to identical
results. Further, magnetic field and temperature are expressed in energy
units, and entropy is expressed in $k_{B}$ units.

To obtain temperature dependences of MCE parameters for the system with
Hamiltonian (1), we used the Metropolis algorithm within the classical Monte
Carlo method. Calculations were carried out for a 64$\times$64 square lattice
with periodic boundary conditions. In our program, we ensured that the
condition for a constant concentration of impurities, i.e.,
$nN=\sum_{i}P_{0i}=const$ , is held, since the variation in the state of an
arbitrarily chosen pair of lattice sites $a$ and $b$ was performed while sum
$P_{0a}+P_{0b}$ was held constant.

We implemented a variant of this algorithm for parallel computations using
NVIDIA graphic accelerators, which enabled us to perform calculations for 64
copies of the system simultaneously. For each system copy, we lowered the
temperature, while keeping the external magnetic field constant, from value
$T_{1}/J=1.0$, which is around twice the magnetic ordering temperature, to
$T_{0}/J=0.01$ at a step of $T/J=0.01$. We then varied the external magnetic
field from $h/J=0$ to $h_{m}/J=0.04$ at a step of $h/J=0.001$. Modeling was
performed for impurity concentration $0<n<0.5$ at a step $n=0.1$. For each
value for $n$, $T$ and $h$ we obtained average values for energy, specific
heat $C$ and magnetization $M$. This allowed us to use the trapezium rule for
discretization of integrals and a threepoint method for derivative $(\partial
M/\partial T)_{h}$ in Eqs.(2)–(4).

With Eqs. (2) and (3) , we can calculate the variation of magnetic entropy
using different Monte Carlo data, i.e., the specific heat and magnetization.
The results of these calculations are in fairly good agree ment in case of
strong coupling ($\tilde{J}>V$) but at variance in case of weak coupling
($\tilde{J}<V$). This is shown in Fig. 2 for impurity concentration $n=0.3$.

Fig. 2:  Magnetic entropy change for $n=0.3$ in the case of (a) strong
coupling ($V/\tilde{J}=0.4$) and (b) weak coupling ($V/\tilde{J}=4$). Results
of calculations carried out using Eqs.(2) and (3) are shown with dotted lines
and filled circles, respectively. A continuous line pertaining to the case of
weak coupling takes into account a contribution calculated by Eq.(6). The
value corresponds to the average concentration of paramagnetic clusters for
$n=0.3$ shown in Fig. 3.

In contrast to the case of strong coupling, with a strong densitydensity
interaction between impurities ($V>\tilde{J}$) , the value of paired
distribution function for the nearest neighbors turns 0 at a sufficiently low
temperature and leads to charge ordering for $n\geq 0.3$. A picture
characteristic of such charge ordering is shown in Fig.3. As a result,
isolated spin clusters surrounded by nonmagnetic impurities are formed. These
clusters behave like paramagnetic centers and furnish additional contribution
$S_{0}$ to the system entropy at the zero temperature. The simplest estimation
of this contribution can be based on the equation

$S_{0}=N_{cl}\ln 2,$ (5)

where $N_{cl}$ is the number of isolated spin clusters. To determine average
number of clusters on a 64$\times$64, lattice, we used a separate program that
generated random distributions of charged impurities having a minimal energy.
The results of calculations are presented in Fig.3. Considering the
contribution given by Eq.(5), formula (2) reads

$\Delta
S_{M}(T,h)=S_{0}+\int\limits_{T_{min}}^{T}\dfrac{C(T^{\prime},h_{m})-C(T^{\prime},0)}{T^{\prime}}dT^{\prime},$
(6)

where $T_{min}$ is the minimum temperature used in Monte Carlo simulations. As
is shown in Fig.2, this improved consistency with $\Delta S_{M}$, M dependence
that was calculated using Eq.(3). Thus, differences in results of magnetic
entropy $\Delta S_{M}$ measurements that were obtained through Eqs.(2) and (3)
may point to some hidden frustrations in the system.

Fig. 3:  (a) Snapshot showing a configuration of the ground state for $n=0.3$
and $V/\tilde{J}=4$ (weak coupling). Lattice sites with $S_{z}=1$ and
$S_{z}=-1$ are represented by white and light-grey squares, respectively.
Positions of nonmagnetic impurities ($S_{z}=0$) are designated with black
squares. (b) Average concentration of isolated paramagnetic clusters on a
64$\times$64 lattice. Inset: the results for clusters comprised of 5 and 8
spins. The difference between maximum values for magnetic entropy variation
$\Delta S_{M}$, as calculated using Eqs.(2) and (3) and expressed in the units
of average concentration of clusters (Eq.(5)), is shown with filled circles.

##  3 Results

The dependence of isolated cluster concentration and the difference between
maximum values for magnetic entropy change $\Delta S_{M}$ as calculated by
Eqs. (2) and (3) and expressed in the units of average concentration, are
presented in Fig.3. We highlight a good agreement between the cluster
concentration and the frustration data for the ground state of the system
obtained from MCE. Partial contribution of clusters consisting of 5 and 8
spins are shown in the inset to this figure. The maximum values for these
contributions are an order of magnitude lower than the total concentration of
isolated paramagnetic clusters for a given $n$, which draws us to a conclusion
that smallest clusters, including 1 spin, make a decisive contribution to the
ground state entropy. In this case, the concentration of smallest paramagnetic
clusters increases monotonically with increasing $n$ reaching maximum
$N_{cl}/N=0.5$ at $n=0.5$.

Fig. 4:  Magnetic entropy change $\Delta S_{M}$ and adiabatic temperature
change $\Delta T_{ad}$ for concentrations $n$ of nonmagnetic impurities in the
range of $00.4$ in case of strong coupling ($V/\tilde{J}=0.4$).

Using a discrete approximation for Eqs. (3) and (4) we calculated isothermal
magnetic entropy change $\Delta S_{M}$ and adiabatic temperature change
$\Delta T_{ad}$ on the basis of Monte Carlo data. The results for $\Delta
S_{M}$ and $\Delta T_{ad}$ in the case of strong coupling, i.e.,
$V/\tilde{J}=0.4$, for $n=0.0,0.1,0.2,0.3,0.4$ are shown in Fig.4. The
temperatures corresponding to maximum values for the two parameters follow
approximately the variation of magnetic ordering temperature with $n$. The
peak width grows, while their height diminishes as $n$ increases, because of
impurity-induced smearing of the phase FM transition. At low temperatures,
dependences $\Delta S_{M}$ and $\Delta T_{ad}$ feature an additional peak
caused by phase separation. $\Delta S_{M}$ and $\Delta T_{ad}$ in the case of
weak coupling, i.e., $V/\tilde{J}=4$, for $n=0.0,0.1,0.2,0.3,0.4,0.5$ are
shown in Fig.5. For $n=0.0$ and $n=0.1$, these results are nearly identical to
those for strong coupling. However, qualitative differences are seen at
$n>0.1$ , especially at low temperatures. As was discussed above, charge
ordering in the impurity subsystem leads to the formation of isolated spin
clusters. This causes a paramagnetic response of the MCE parameters that gives
maximum value $\Delta T_{ad}/J\approx 0.55$ for $n=0.5$ and $T/J=0.14$.

Fig. 5:  Magnetic entropy change $\Delta S_{M}$ and adiabatic entropy change
$\Delta T_{ad}$ for concentrations $n$ of nonmagnetic impurities in the range
of $00.5$ in case of weak coupling ($V/\tilde{J}=4$). Inset: results for
$n=0.5$.

##  4 Conclusions

We performed Monte Carlo simulations of a twodimensional Ising system
containing a fixed concentration of mobile nonmagnetic charged impurities,
which enabled us to construct temperature dependences of the MCE parameters.
It was shown that in case of weak coupling, for nonmagnetic impurity
concentration $n>0.1$ , isolated spin clusters are formed in the ground state
of the system. This leads to nonzero entropy of the systems ground state. The
feasibility to detect a frustration in the ground state is demonstrated using
the magnetic entropy variation data. The MCE parameters were calculated for
the cases of strong and weak coupling.

## Funding

The work was supported within the Competitiveness Enhancement Program of the
Ural Federal University (Russian Federation Government Act 211, agreement no.
02.A03.21.0006), by the Ministry of Education and Science, Russian Federation
(project FEUZ$-2020-0054$), and by the Russian Foundation for Basic Research
(research project no. $18-32-00837\backslash 18$).

## Conflict of interest

The authors declare that they have no conflicts of interest.

## References

  * [1] M. E. Zhitomirsky. Phys. Rev. B 67, 104421 (2003).
  * [2] M. E. Zhitomirsky, A. Honecker. J. Stat. Mech.: Theory Exp. 2004, P07012 (2004).
  * [3] A. Honecker, S. Wessel. Physica B 378-380, 1098 (2006).
  * [4] B. Schmidt, P. Thalmeier, N. Shannon. Phys. Rev. B 76, 125113 (2007).
  * [5] J. S. Amaral, J. N. Goncalves, V. S. Amaral. IEEE Trans. Magn. 50, 1 (2014).
  * [6] M. Žukovič. J. Magn. Magn. Mat. 374, 22 (2015).
  * [7] S. Katsura, B. Tsujiyama. Ferro- and Antiferromagnetism of Dilute Ising Model, in: C. Domb (Ed.), Proceedings of the Conference on Phenomena in the Neighborhood of Critical Points, National Bureau of Standards, Washington, D.C., 1965, pp. 219–224.
  * [8] M. Blume, V. J. Emery, R. B. Griffiths. Phys. Rev. A 4, 1071 (1971).
  * [9] A. S. Moskvin. J. Phys.: Condens. Matter 25, 085601 (2013).
  * [10] B. G. Shen, J. R. Sun, F. X. Hu, H. W. Zhang, Z. H. Cheng. Adv. Mater. 21, 4545 (2009).
  * [11] A. A. Inishev, E. G. Gerasimov, N. V. Mushnikov, P. B. Terent’ev, V. S. Gaviko. Physics of Metals and Metallography 119, 1036 (2018).
  * [12] H. Zhang, R. Gimaev, B. Kovalev, K. Kamilov, V. Zverev, A. Tishin. Physica B 558, 65 (2019).
  * [13] Y. D. Panov, A. S. Moskvin, A. A. Chikov, I. L. Avvakumov. J. Low Temp. Phys. 185, 409 (2016).
  * [14] Y. D. Panov, A. S. Moskvin, A. A. Chikov, I. L. Avvakumov. Journal of Superconductivity and Novel Magnetism 29, 1077 (2016).
  * [15] Y. D. Panov, K. S. Budrin, A. A. Chikov, A. S. Moskvin. JETP Letters 106, 440 (2017).
  * [16] Y. Panov, V. Ulitko, K. Budrin, A. Chikov, A. Moskvin. J. Magn. Magn. Mat. 477, 162 (2019).