The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
33,507 characters
Stochastic Frontier meets Breakdown Frontier
\nonstopmode
\title[Stochastic Frontier meets Breakdown Frontier]{Stochastic Frontier meets Breakdown Frontier}
\author{Santiago Acerenza, Francisco Rosas}
\address{Universidad ORT Uruguay}
\noindent \date{\scriptsize{ The present version is as of \today.
}}
\begin{abstract}
This paper studies sensitivity analysis of Stochastic Frontier Models. We elaborate relaxations of the baseline assumptions in the Stochastic Frontier Models and characterize the identified set under this relaxations. Furthermore, we derive the breakdown frontier for a relevant parameter of interest, the average inefficiency of a production unit. We show an application of the procedures on a well known dataset, and make the code available for the interested practitioner.
\end{abstract}
\maketitle
\maketitle
{\footnotesize \textbf{Keywords}: Stochastic Frontier Analysis, Technical inefficiency, Breakdown Frontier, Partial Identification.
\textbf{JEL subject classification}: C18, C2, C23, C51, D24 }
\section{Introduction}
In this paper, we aim at responding the question, what are the maximum relaxations of the imposed assumptions on the unobservables which still sustain our conclusion in terms of efficiency of production units, given the observed data? Since its first introduction by \cite{aigner1977formulation} and \cite{meeusen1977efficiency}, stochastic frontier analysis has been applied to study the productivity and efficiency of production units in various economic sectors, such as banking \citep{ferrier1990measuring,adams1999semiparametric,kumbhakar2005measuring,malikov2016cost}, healthcare \citep{zuckerman1994measuring,rosko2001cost,greene2004distinguishing,mutter2013investigating,comans2020cost} and agriculture \citep{Bravo2007,Bravo2017, Trestini, Qushim2013,Ozden,Otieno,Gatti,Nwigwe, Qushim2018,Martinez,xin2025environmental,Lanfranco2013,Garcia2019,Garcia2022,Aguirre2024a,Aguirre2024b}. The use of stochastic frontier models is still present in the current empirical practice, and significant extensions have been made to the methodology since the seminal works from \cite{aigner1977formulation} and \cite{meeusen1977efficiency} to account for more specific contexts, exploit more information, and to obtain more robust results. To name a few,\cite{schmidt1984production,greene2005fixed, colombi2014closed} and \cite{kumbhakar2014technical} contributed in the realm of panel data model extensions, including the estimation of the persistent and transient components of productive inefficiency. \cite{cornwell1990production,wang2002one,caudill1995frontier} extended the analysis for models with determinants of inefficiency. \cite{simar2017nonparametric,wang2024flexible,centorrino2024nonparametric,zheng2024robust} expanded the methodology for semi-parametric and non-parametric stochastic frontier models (SFM). See \cite{nguyen2022efficiency} for a thorough review as well as software implementation details.
\par
Although these significant advances had made conclusions related to (in)efficiency become more robust to miss-specification, no sensitivity analysis per-se has been introduced in the literature. For sensitivity analysis \citep{masten2020inference} we refer to a multidimensional approach that evaluates the robustness of empirical conclusions to simultaneous relaxations of a set of identifying assumptions. Typical approaches start with assumptions to identify a parameter, however, this approach begins with a specific conclusion and determines the weakest set of assumptions required for that conclusion to hold given the observed data. This is formalized in a breakdown frontier, which maps the tradeoffs between relaxing different assumptions while still maintaining the validity of the desired conclusion.
\par
This is relevant, firstly, because all the previous mentioned papers require certain identifying assumptions (no matter how slack they are) to provide a conclusion about (in)efficiency. Secondly, this analysis determines the threshold at which a straightforward model remains robust, tolerating sufficient deviation from baseline assumptions. Consequently, researchers may retain simpler models, maximizing both interpretability and implementation efficiency.
\par
Following \cite{masten2020inference} we bridge this gap introducing a Breakdown frontier analysis in this context. The breakdown frontier generalizes the concept of an “identification breakdown point” introduced by \cite{horowitz1995identification}, a one dimensional breakdown frontier. As noted by \cite{masten2020inference} the breakdown frontier approach requires six main steps: (a) specifying a parameter of interest (typically in our context the inefficiency for a particular production unit), (b) specifying a set of baseline assumptions (for example truncated normality of the inefficiency term), (c) defining a class of assumptions indexed by a sensitivity parameter, which delivers a nested sequence of identified sets with the baseline assumptions obtained at one extreme and the no assumptions bounds obtained at the other, (d) characterizing identified sets for the parameter of interest as a function of the sensitivity parameter, (e) using those identified sets to define the breakdown frontier for a conclusion of interest, and (f) developing estimation and inference procedures for that frontier based on its characterization.
\par
We will follow these steps with particular focus on the assumptions from the seminal model from \cite{aigner1977formulation} but this approach can be extended with proper modifications to the other models mentioned in this introduction above.
\section{Framework}
\subsection{Basic Set-up and the parameter of interest}
Let $y_i$ be the logarithm of an output of interest. Then, the classic model for stochastic frontiers states:
\begin{eqnarray*}
y_i=f(\theta_x,x_i)+v_i-e_i
\end{eqnarray*}
Where $f(\theta_x,x_i)$ is a known function of the inputs $x_i$ up to parameters $\theta_x \in \mathbb{R}^p$, $v_i$ is an error term, and $e_i$ is the unobservable component associated with technical (in)efficiency, where $e_i\geq 0$. It is common, in the stochastic frontier literature, to try to identify:
\begin{eqnarray*}
E[exp^{-e_i}|v_i-e_i]
\end{eqnarray*}
If one is interested in the estimates of individual (in)efficiency of a production unit, this measure captures it.
\subsection{The structure of distributional assumptions and its relaxations}
It is also common in the literature to impose additional structure on $v_i$ and $e_i$ to identify this parameter of interest. A typical set of assumptions is normality of $v_i$ and some form of truncated distribution on $e_i$, potentially with mean restrictions and heterosckedasticity.
\par
The intuition behind the sensitivity analysis is the following. Assume that there exist two unknown constants $b$ and $c$ that, respectively, measure the distance between the assumed distribution of $v$ and $e$ to the underlying or true distribution of these error terms. Then, in the context of the model presented above, suppose that the average value of technical inefficiency for a particular $u$ was estimated at $0.5$. The objective of the sensitivity analysis is to find the combinations of $c$ and $b$, which actually imply violations of the assumptions, such that $0.5$ is still a robust result. This is the so-called robust region (RR), in other words, a measure of how far the true distribution has to be from the imposed assumption in order to invalidate the estimated conclusion. The objective of the sensitivity analysis is to find, conditional on the estimated $0.5$ and a given value of $c$ that can be conceived as a tolerable violation of the efficiency measure, which is the value of the unknown constant $b$. This is the breakdown frontier, that is graphically depicted in Figure \ref{fig1}.
\begin{figure}[H]
\includegraphics[scale=0.5]{convex_breakdown_frontier.pdf}
\caption{The Breakdown Frontier. Source: Own elaboration based on \cite{masten2020inference}}
\label{fig1}
\end{figure}
The blue region are the set of combinations of $c$ and $b$, i.e. violations of the distributional assumptions, for which the conclusion of average (in)efficiency equal to some level holds. The yellow regions are the combinations for which we have not sufficient evidence in favor of the conclusions, given the imposed assumptions. The negative slope of the frontier implies that there is a tradeoff in reducing the violation of one of the assumptions, either on the error term ($v$) or in the efficiency term ($e$) that implies a larger violation of the other assumption. Or what is the same, if a relaxation of one of the assumptions is permitted, this comes at the benefit of a lower violation of the other assumption. The nonlinearity of the decreasing curve means that for a given relaxation of one assumption, say on $b$, it is required a reduction in the violation of the other ($c$) for the conclusions to hold, that varies depending on the value of $b$. Its convexity implies that a relaxation of $b$ when we allowed a small violation of that assumption permits a sizable reduction in the violation of the other assumption. And conversely, the same relaxation of $b$ when it is already broadly relaxed only allows for a small gain in the reduction of the violation of $c$.
Our starting point is to assume a relaxation of these distributional assumptions as a function of a parameter that measures the distance between these typically imposed baseline assumptions and the no-assumptions (and no identification of the parameter of interest) case.
\begin{assumption}[Parametric distance]\label{AsPara}
Let $f_{e_i}(e)$ denote the density function of the variable $e_i$ and let $f_{v_i|e_i}(v|e)$ denote the conditional distribution of $v_i$ given $e_i$. Furthermore, let $f_{e_i}(e;\theta_e), f_{v_i|e_i}(v|e; \theta_v)$ denote parametric density functions specified by the researcher. Then:
\begin{eqnarray*}
\sup_{e} |f_{e_i}(e)-f_{e_i}(e;\theta_e)|&\leq& c \\
\sup_{v} |f_{v_i|e_i}(v|e)-f_{v_i|e_i}(v|e; \theta_v)|&\leq& b \quad \forall e\\
\end{eqnarray*}
\end{assumption}
Note that the case that $c=0,b=0$ would be the classic stochastic frontier model, which in the case that $v_i$ follows a Normal distribution and $e_i$ a truncated normal, we are under the case of \cite{aigner1977formulation}. Moving away from $c=0$ and $b=0$ implies relaxing the structure of the distributional assumptions of the model.
\subsection{The identified set for the parameter of interest}
In this section, we present how the identified set for the parameter of interest looks like under the relaxation of the classic assumptions ($b=0,c=0$). This result is collected in the following lemma.
\begin{lemma}\label{LemmaIdentifiedset}
Under Assumption \ref{AsPara}:
\begin{eqnarray}\label{EqBP5}
E[exp^{-e_i}|u] &\leq& \frac{ \int_{0}^{\infty} exp^{-e}f_{v_i|e_i}(u+e|e;\theta_v) f_{e_i}(e;\theta_e)de+ c\int_{0}^{\infty} exp^{-e}f_{v_i|e_i}(u+e|e;\theta_v)de}{\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)f_{e_i}(e; \theta_e)de- c\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)de-b} \nonumber \\
&+& \frac{b \int_{0}^{\infty} exp^{-e}f_{e_i}(e;\theta_e)de+bc}{\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)f_{e_i}(e; \theta_e)de- c\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)de-b} \nonumber \\
E[exp^{-e_i}|u] &\geq& \frac{ \int_{0}^{\infty} exp^{-e}f_{v_i|e_i}(u+e|e;\theta_v) f_{e_i}(e;\theta_e)de- c\int_{0}^{\infty} exp^{-e}f_{v_i|e_i}(u+e|e;\theta_v)de}{\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)f_{e_i}(e; \theta_e)de+ c\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)de+b} \nonumber \\
&+& \frac{ -b \int_{0}^{\infty} exp^{-e}f_{e_i}(e;\theta_e)de+bc}{\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)f_{e_i}(e; \theta_e)de+ c\int_{0}^{\infty} f_{v_i|e_i}(u+e|e;\theta_v)de+b}
\end{eqnarray}
\end{lemma}
Where $u=v-e$.
\par
Suppose that the average value of inefficiency for a particular $u$ from the model with $c,b=0$ was $0.5$, then this can be used to find the combinations of $c,b$ (of violations of the assumptions) such that $0.5$ is a robust result. This is the robust region, which would be a measure of how far away has to be the true values from the imposed assumption in order to reject our observed conclusion. This yields a measure of sensitivity.
\subsection{An application to the Normal-Truncated Normal Model: Using the identified set to define the breakdown frontier for a conclusion of
interest}
Suppose that $e_i,v_i$ are independent of each other and also assume that $e_i\sim TN(\mu,\sigma^2_e)$ and $v_i\sim N(0,\sigma^2_v)$. Then:
\begin{eqnarray*}
f_{e_i}(e,\theta_e)&=& \frac{1}{\sqrt{2\pi\sigma_e^2}}\frac{exp^{\frac{-(e-\mu)^2}{2 \sigma^2_e}}}{1-\Phi(\frac{-\mu}{\sigma_e})} \\
f_{v_i|e_i}(v|e;\theta_v)&=& \frac{1}{\sqrt{2\pi\sigma_v^2}} exp^{\frac{-(v)^2}{2 \sigma^2_v}}
\end{eqnarray*}
Define then for the previously mentioned densities,
$$ \int_{0}^{\infty } exp^{-e} f_{v_i|e_i}(u+e|e;\theta_v)de \equiv \Delta_1(u;\sigma_v),
\int_{0}^{\infty } exp^{-e} f_{v_i|e_i}(u+e|e;\theta_v)f_{e_i}(e,\theta_e)de \equiv \Delta_1(u;\sigma_v,\sigma_e,\mu),$$
$$\int_{0}^{\infty } exp^{-e}f_{e_i}(e,\theta_e)de \equiv \Delta_1(u;\sigma_e,\mu),
\int_{0}^{\infty } f_{v_i|e_i}(u+e|e;\theta_v)de \equiv \Delta_2(u; \sigma_v),$$
$$\int_{0}^{\infty } f_{v_i|e_i}(u+e|e;\theta_v)f_{e_i}(e,\theta_e)de \equiv \Delta_2(u;\sigma_v,\sigma_e,\mu).$$
Where their closed form expressions can be found as a result of Lemma \ref{LemmaAuxiliar}.
\subsubsection{The identified set}
By an application of equation \ref{EqBP5}:
\begin{eqnarray}\label{NormalIdentifiedset}
E[exp^{-e_i}|v_i-e_i=u] &\leq& \frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)+c \Delta_1(u; \sigma_v)+b\Delta_1(u;\sigma_e, \mu)+bc}{\Delta_2(u; \sigma_v,\sigma_e, \mu)-c \Delta_2(u; \sigma_v)-b} \nonumber \\
E[exp^{-e_i}|v_i-e_i=u] &\geq& \frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)-c \Delta_1(u; \sigma_v)-b\Delta_1(u;\sigma_e, \mu)+bc}{\Delta_2(u; \sigma_v,\sigma_e, \mu)+c \Delta_2(u; \sigma_v)+b}
\end{eqnarray}
\subsubsection{The breakdown frontier}
We now formally define the breakdown frontier in this context. We also introduce the robust region, the area below the breakdown frontier.
\par
Suppose we run the usual Stochastic Frontier Model (with $c=0,b=0$) and we obtain that $E[exp^{-e_i}|v_i-e_i=u]=e_o$, we then for the breakdown frontier analysis begin with the conclusion that $E[exp^{-e_i}|v_i-e_i=u]\geq e_o$ for a fixed $e_o \in [0,1]$ (since the inefficiency is between 0 and 1). Relative to the baseline assumptions, what are the weakest assumptions that allow us to draw this conclusion, given the observed distribution of the data? Specifically, since larger values of $c$ and $b$ correspond to weaker assumptions, what are the largest values of $c$ and $b$ such that we can still definitively conclude that $E[exp^{-e_i}|v_i-e_i=u]\geq e_o$? This can be answered in two steps following \cite{masten2020inference}. First, it requires gathering all values of $c$ and $b$ such that the conclusion holds. This is called the robust region. Since the lower bound of the identified set by Lemma \ref{LemmaIdentifiedset} and Equation \ref{NormalIdentifiedset} is:
$$ E[exp^{-e_i}|v_i-e_i=u] \geq \frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)-c \Delta_1(u; \sigma_v)-b\Delta_1(u;\sigma_e, \mu)+bc}{\Delta_2(u; \sigma_v,\sigma_e, \mu)+c \Delta_2(u; \sigma_v)+b} $$
And since $c$ and $b$ can take any number between $0$ and plus infinity, we have that in our context the robust region is:
\begin{eqnarray}\label{RobustRegion}
RR=\Big\{c,b \in [0,\infty]^2: \frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)-c \Delta_1(u; \sigma_v)-b\Delta_1(u;\sigma_e, \mu)+bc}{\Delta_2(u; \sigma_v,\sigma_e, \mu)+c \Delta_2(u; \sigma_v)+b}\geq e_0 \Big\}
\end{eqnarray}
Note that when $c=0,b=0$, $e_o=\frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)}{\Delta_2(u; \sigma_v,\sigma_e, \mu)}$ and since the previous function of $c,b$ is increasing in both of its arguments $c=0,b=0$ is included in the robust region and thus is not empty.
\par
The breakdown frontier is the set of points on the boundary of the robust region. Namely:
\begin{eqnarray}\label{Breakdownfrontier}
BF(e_o)=\Big \{ c,b \in [0,\infty]^2: \frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)-c \Delta_1(u; \sigma_v)-b\Delta_1(u;\sigma_e, \mu)+bc}{\Delta_2(u; \sigma_v,\sigma_e, \mu)+c \Delta_2(u; \sigma_v)+b}=e_0\Big\}
\end{eqnarray}
In order to express this in a simple graphical manner, we would like to see what are the values of $c$ for any value of $b$ that is in the frontier. This would be an equivalent way of characterizing the previous set. Thus, since the value of $e_o$ is fixed, solving for $b$:
\begin{eqnarray}\label{bcenodu}
b(c,e_o)=\frac{c[\Delta_1(u; \sigma_v)+e_o\Delta_2(u; \sigma_v)]+e_o\Delta_2(u; \sigma_v,\sigma_e, \mu)-\Delta_1(u; \sigma_v,\sigma_e, \mu)}{c-\Delta_1(u;\sigma_e, \mu)-e_o}
\end{eqnarray}
This allows us to get the following analytical expression for the $BF(e_o)$:
\begin{eqnarray}
BF(e_o)= \max\{ b(c,e_o),0\}
\end{eqnarray}
\subsection{The Normal-Truncated Normal Model continued: Estimation and Inference procedures for that frontier}
In this section we study estimation and inference on the breakdown frontier defined above. The breakdown frontier is a known functional of the $\Delta$ objects which are known or estimable using classic maximum likelihood estimation as it is the case in \cite{aigner1977formulation} model. Hence we propose plug-in estimators of the breakdown frontier. We can get $\sqrt{n}$-consistency and asymptotic distributional results using population smoothing procedures\footnote{ We can alternatively use the delta method for directionally differentiable functionals \cite{fang2019inference} and implement a bootstrap procedure to construct asymptotically valid lower confidence bands for the breakdown frontier following the numerical delta method from \cite{hong2018numerical}. \cite{Masten2019} denotes that there is no clear dominance between the numerical delta method and the population smoothing.} that recovers an approximation of this maximum that is not only directionally differentiable but fully differentiable, and therefore, standard asymptotic normality is obtained, and thus, the classic Nonparametric Bootstrap can be performed \citep{fang2019inference}.
\par
Note that the main elements of the sensitivity frontier are the elements that can be estimated fron the "Null" model, a.k.a the model with $c=0,b=0$. Thus, we need to estimate $\theta=(\theta_x,\mu,\sigma_v,\sigma_e) \in \mathbb{R}^{p+3}$. This can be done via optimizing the following log-likelihood (for an $i.i.d$ sample of size $n$):
\begin{eqnarray}\label{LogLikelihood}
LnL&=&\sum_{i=1}^n ln f(y_i,x_i;\theta)=\sum_{i=1}^n \Bigg( -\frac{1}{2}ln(2\pi)-ln\Big( \sqrt{(\sigma_v^2+\sigma_e^2)}\Big)-ln\Phi\Big( \frac{\mu}{\sigma_e} \Big) \nonumber \\
&+& ln\Phi\Big( \frac{(1-\frac{\sigma^2_e}{\sigma^2_e+\sigma^2_v}
)\mu- \frac{\sigma^2_e}{\sigma^2_e+\sigma^2_v}
(y_i-f(\theta_x,x_i))}{\sqrt{(\sigma^2_e+\sigma^2_v)\frac{\sigma^2_e}{\sigma^2_e+\sigma^2_v}(1-\frac{\sigma^2_e}{\sigma^2_e+\sigma^2_v})}} \Big)-\frac{1}{2}\Big(\frac{y_i-f(\theta_x,x_i)+\mu}{ \sqrt{(\sigma_v^2+\sigma_e^2)}} \Big)^2 \Bigg)
\end{eqnarray}
Where we obtain estimates for $\theta=(\theta_x,\mu,\sigma_v,\sigma_e)$ namely $\widehat{\theta}_n=(\widehat{\theta}_{xn},\widehat{\mu}_n,\widehat{\sigma}_{vn},\widehat{\sigma}_{en})$. Note that $ f(y_i,x_i;\theta)$ is the likelihood of observation $i$.
By standard asymptotic results,
\begin{eqnarray*}
\sqrt{n}\begin{pmatrix}
\theta-\widehat{\theta}_n
\end{pmatrix}&\xrightarrow{d} N(\boldsymbol{0}, \mathcal{H}_{\theta})&
\end{eqnarray*}
Where $ \mathcal{H}_{\theta}=-E[\frac{\partial ^2}{\partial \theta \partial \theta'} lnf(y_i,x_i;\theta)]$ or in more detail,
\begin{eqnarray*}
\mathcal{H}_{\theta}=-E \begin{pmatrix}
\frac{\partial ^2}{\partial \theta_{x1} \partial \theta_{x1}} lnf(y_i,x_i;\theta)
& \cdots & \frac{\partial ^2}{\partial \theta_{x1} \partial \theta_{xp}} lnf(y_i,x_i;\theta) & \cdots & \frac{\partial ^2}{\partial \theta_{x1} \partial \sigma_{e}} lnf(y_i,x_i;\theta) \\\
\vdots & \vdots& \vdots & \vdots & \vdots \\
\frac{\partial ^2}{\partial \sigma_{e} \partial \theta_{x1} } lnf(y_i,x_i;\theta) & \cdots & \frac{\partial ^2}{ \partial \sigma_{e} \partial \theta_{xp}} lnf(y_i,x_i;\theta) & \cdots & \frac{\partial ^2}{ \partial \sigma_{e} \partial \sigma_{e}} lnf(y_i,x_i;\theta)
\end{pmatrix} \quad
\end{eqnarray*}
The variance covariance can be estimated consistently with $ \widehat{\mathcal{H}}_{\theta n}=-\frac{1}{n}\sum_{i=1}^n [\frac{\partial ^2}{\partial \widehat{\theta}_n \partial\widehat{\theta}_n'} lnf(y_i,x_i;\widehat{\theta}_n)]$.
\par
Note that the $\Delta$´s are all continuously differentiable functions of $\theta$, thus all of them have an asymptotically Gaussian behavior. Concretely, if
$$\Delta=(\Delta_1(u;\sigma_v),\Delta_2(u;\sigma_v), \Delta_2(u;\sigma_v,\sigma_e,\mu), \Delta_1(u;\sigma_v,\sigma_e,\mu), \Delta_1(u;\sigma_e,\mu)) \in \mathbb{R}_+^{5}$$
And
$$\widehat{\Delta}_n=(\Delta_1(u;\widehat{\sigma}_{vn}),\Delta_2(u;\widehat{\sigma}_{vn}), \Delta_2(u;\widehat{\sigma}_{vn},\widehat{\sigma}_{en},\widehat{\mu}_n), \Delta_1(u;\widehat{\sigma}_{vn},\widehat{\sigma}_{en},\widehat{\mu}_n), \Delta_1(u;\widehat{\sigma}_{en},\widehat{\mu}_n))$$
Then,
\begin{eqnarray*}
\sqrt{n}(\Delta-\widehat{\Delta}_n)\xrightarrow{d} N(\boldsymbol{0},\boldsymbol{D}_\Delta \mathcal{H}_{\theta} \boldsymbol{D}_\Delta')
\end{eqnarray*}
Where,
\begin{eqnarray*}
\boldsymbol{D}_\Delta = \begin{pmatrix}
\frac{\partial}{\partial \theta_{x1}} \Delta_1(u;\sigma_v) & \cdots & \frac{\partial}{\partial \sigma_e} \Delta_1(u;\sigma_v) \\
\vdots & \vdots & \vdots \\
\frac{\partial}{\partial \theta_{x1}} \Delta_1(u;\sigma_e,\mu) & \cdots & \frac{\partial}{\partial \sigma_e} \Delta_1(u;\sigma_e,\mu)
\end{pmatrix}
\end{eqnarray*}
Furthermore, as long as $c-\Delta_1(u;\sigma_e,\mu)-e_o \neq 0$ it will then, by another application of the delta method, also be true that:
\begin{eqnarray}\label{EqAsymptoticb}
\sqrt{n}(b(c,e_o)-\widehat{b}(c,e_o)_n) \xrightarrow{d} N(0, \sigma_{b(c,e_o)}^2)
\end{eqnarray}
Where,
\begin{eqnarray*}
\sigma_{b(c,e_o)}^2&=& \boldsymbol{D}_{b(c,e_o)} \boldsymbol{D}_\Delta \mathcal{H}_{\theta} \boldsymbol{D}_\Delta' \boldsymbol{D}_{b(c,e_o)} ' \\
\boldsymbol{D}_{b(c,e_o)} &=& \begin{pmatrix}
\frac{\partial}{\partial \Delta_1(u;\sigma_v) } b(c,e_o) & \cdots & \frac{\partial}{\partial \Delta_1(u;\sigma_e,\mu) } b(c,e_o)
\end{pmatrix}
\end{eqnarray*}
The next step is to obtain an asymptotic distribution for an estimator of $BF(e_o)$, but this is a non-differentiable function and thus the standard results break down. Here is where population smoothing enters.
\par
Population smoothing consists of using the smooth approximation of the maximum, that will allow to obtain a differentiable approximation to which we can apply standard bootstrap results. This approach is advantageous because we can use a non-parametric bootstrap to construct confidence intervals for the breakdown frontier. Consequently, we avoid estimating complex, convoluted variance components and bypass the need to handle a non-pivotal statistics.
\begin{eqnarray*}
\max\{m_1,..m_K\} \geq \frac{\sum_{k=1}^{K}m_{k}e^{\rho m_{k}}}{\sum_{k=1}^{K}e^{\rho m_{k}}}
\end{eqnarray*}
Where as $\rho \rightarrow \infty$, $\frac{\sum_{k=1}^{K}m_{k}e^{\rho m_{k}}}{\sum_{k=1}^{K}e^{\rho m_{k}}} \rightarrow \max\{m_1,..m_K\}$.
\par
Thus, in our context, we focus on:
\begin{eqnarray*}
\text{soft}BF(e_o)=\frac{b(c,e_o) e^{\rho b(c,e_o)}}{1+e^{\rho b(c,e_o)}}
\end{eqnarray*}
This is a differentiable function and thus, asymptotically normal point-wise in $c$ (or uniform under extra regularity conditions), then confidence intervals can be constructed by relying on the now asymptotic normality of:
\begin{eqnarray}\label{Dist_softBF}
\sqrt{n}\Big(\frac{b(c,e_o) e^{\rho b(c,e_o)}}{1+e^{\rho b(c,e_o)}} -\frac{\widehat{b}(c,e_o)_n e^{\rho \widehat{b}(c,e_o)_n}}{1+e^{\rho \widehat{b}(c,e_o)_n}} \Big) \xrightarrow{d} N(0, \sigma^2_{\text{soft}BF(e_o)})
\end{eqnarray}
\begin{eqnarray*}
\sigma^2_{\text{soft}BF(e_o)}=\Big( \frac{\partial}{\partial b(c,e_o)} \frac{b(c,e_o) e^{\rho b(c,e_o)}}{1+e^{\rho b(c,e_o)}} \Big)^2 \sigma_{b(c,e_o)}^2
\end{eqnarray*}
In appendix \ref{AppInfluencefunctions} we collect the influence function of our estimator which can then be used if of interest to construct a de-biased estimator of the form of:\footnote{See \cite{hines2022demystifying,fisher2021visually}.}
$$ \frac{\widehat{b}(c,e_o)_n e^{\rho \widehat{b}(c,e_o)_n}}{1+e^{\rho \widehat{b}(c,e_o)_n}} +\frac{1}{n}\sum_{i}^n \widehat{\psi}_{\frac{b(c,e_o) e^{\rho b(c,e_o)}}{1+e^{\rho b(c,e_o)}}}(x_i,y_i).$$
Where $\widehat{\psi}$ is the estimator of the influence function of the estimator.
\subsubsection{What happens when $\rho$ grows?}
If we pick $\rho$ to be a function of the data that grows at a rate slower than the rate our estimators converge, we can in a similar spirit as \cite{Bontempsetal2024} still obtain a consistent statistic for the maximum\footnote{ In fact we can construct a pivotal test statistic to invert and obtain confidence regions this is a result from \cite{Bontempsetal2024}. Up to \today \quad the results from \cite{Bontempsetal2024} are not publicly available but we had knowledge of them from a presentation in IAAE-2024.} Thus, we can claim:
\begin{lemma}\label{LemmaSmooth}
As long as $\rho_n\rightarrow \infty$ grows slower than $(\widehat{b}(c,e_o)_n-b(c,e_o)) \rightarrow 0$ then:
\begin{eqnarray*}
\Big(\max\{0,b(c,e_o)\} -\frac{\widehat{b}(c,e_o)_n e^{\rho_n \widehat{b}(c,e_o)_n}}{1+e^{\rho_n \widehat{b}(c,e_o)_n}} \Big)\xrightarrow{p} 0
\end{eqnarray*}
\end{lemma}
\section{Empirical Illustration}
In this section, we provide a empirical illustration of the tool proposed in the previous sections, following the data used in \cite{nguyen2022efficiency}, consisting of output and input information about rice producers in the Philippines.\footnote{Downloaded from the built in package \url{https://frontier.r-forge.r-project.org/}. Details of the data can be found in \url{https://www.rdocumentation.org/packages/frontier/versions/1.1-8/topics/riceProdPhil}.}
The dataset includes the information about 43 rice producers in Tarlac, Philippines from 1990 to 1997. Out of the original dataset, we use one output $y_i$ (freshly threshed rice in tones) and three inputs $x_i$, the area planted (in hectares), labor used (man-days of family and hired labor), and fertilizer used (active ingredients in kilograms).
We run the SPF model with the output specification as $y_i=f(\theta_x,x_i)+v_i-e_i$, where the error term $v_i$ is assumed $v_i\sim N(0,\sigma^2_v)$, and the technical (in)efficiency $e_i$ is $e_i\sim TN(\mu,\sigma^2_e)$. It is also assumed that they are independent of each other. With the estimated frontier and for an arbitrary value of the composite error term $u_i=v_i-e_i$, say $u_i=0$, i.e. at the sample average, we compute the estimated efficiency as $E[exp^{-e_i}|v_i-e_i]=e_0=\frac{\Delta_1(u; \sigma_v,\sigma_e, \mu)}{\Delta_2(u; \sigma_v,\sigma_e, \mu)}=0.8872$. Then, allowing $b$ and $c$ to vary between zero and plus infinity, we can compute the breakdown frontier following equation \ref{Breakdownfrontier}, which uses the estimated values of our $\sqrt{n}-$consistent plug-in estimators $\widehat{\Delta}_n$. Or equivalently, using equation \ref{bcenodu} and given $e_0=0.8872$, we find the values of $b(c,e_0)$ for each proposed value in the range of relaxations of $c$. The resulting breakdown frontier is presented in Figure \ref{fig2}.
\begin{figure}[H]
\includegraphics[scale=0.5]{Breakdown_Frontier.png}
\caption{Breakdown Frontier for $E(exp^{-e}|u=0)$}
\label{fig2}
\end{figure}
The downward slopping curve shows the tradeoffs between the relaxations of the distributional assumption of the efficiency $e$ and of the error term $v$, that maintain the conclusions ($e_0=0.8872$) to hold. The area under the curve is the robust region given in equation \ref{RobustRegion}, which is interpreted as all the combinations of violations of the imposed assumptions for which the mentioned conclusion hold. It is of particular interest to note that violations of the assumption on $e_i$ refers to a miss-specification of the distribution of the true (in)efficiency, however, violations on the assumption of $v_i$ includes both a miss-specification of the distribution but also a violation of the independence between the inefficiency and the error term.
The second part of the discussion in section 2.5 boils down to constructing confidence intervals for the breakdown frontier, allowed by the population smoothing procedures and performed with a nonparametric bootstrapping technique. In particular, we obtain confidence intervals for the breakdown frontier characterized by equation \ref{bcenodu}, which required the use of population smoothing to sort the non-differentiability of this function. The confidence intervals are then constructed relying on the asymptotically normal distribution of $softBF(e_0)$ given in equation \ref{Dist_softBF}, and by means of bootstrapping.
\begin{figure}[H]
\includegraphics[scale=0.5]{Breakdown_FrontierCI.png}
\caption{Smooth Breakdown Frontier for $E(exp^{-e}|u=0)$, pointwise 95 percent confidence intervals. 100 bootstrap replications and $\rho=10$.}
\label{fig3}
\end{figure}
\section{Further remarks}
\subsection{Allowing for heterosckedasticity}
Heterosckedasticity is easily incorporated since the null model can be estimated with different assumptions on the variances and the mean of the errors.
\subsection{Relaxing Exogeneity of the inputs}
The previous analysis focused on the distributional assumptions, but the breakdown frontier can be extended to incorporate relaxations of input exogeneity. A way of doing this would be to incorporate bounds on the following distances $|E(e_i|x)-E(e_i)|, |V(e_i|x)-V(e_i)|$.
\subsection{Non analytical solutions to the main components of the bounds}
The analysis has been exemplified in the case $\Delta$ can be analytically derived, but there will be cases in which expressions like $\int_{0}^{\infty} exp^{-e}f_{v_i|e_i}(u+e|e;\theta_v) f_{e_i}(e;\theta_e)de$ will not have a closed form solution (or the closed form solutions are not easy to implement). In such cases, we could leverage approximation results and obtain estimates for $\Delta$. For example if we obtain an estimate for $\theta_v,\theta_e$, $\widehat{\theta}_v,\widehat{\theta}_e$ we could generate via simulations a sample of size $S$ of realizations from $f_{v_i,e_i}(v,e;\widehat{\theta}_v, \widehat{\theta}_e)=f_{v_i|e_i}(v|e;\widehat{\theta}_v)f_{e_i}(e;\widehat{\theta}_e)$ and estimate $\int_{0}^{\infty} exp^{-e}f_{v_i|e_i}(u+e|e;\theta_v) f_{e_i}(e;\theta_e)de$ with $\sum_{s=1}^{S}exp^{-e_s}f_{v_i|e_i}(u+e_s|e_s;\widehat{\theta}_v) f_{e_i}(e_s;\widehat{\theta}_e)$. Under extra regularity conditions on the simulated sample we could (properly adjusting the variances) obtain similar asymptotic results.
\subsection{Extension to Panel data models}
The empirical application focuses mainly in the cross-sectional setting. Nevertheless, our procedure can be easily adapted to the panel data case in two directions. Firstly, the cross sectional analysis can be implemented year by year and thus construct a set of breakdown frontiers for each year. Secondly, in the presence of fixed effects, where the equation becomes $ y_{i,t}=f(\theta_x,x_{i,t})+v_{i,t}-e_{i,t}+c_i$ our model can be implemented on first differences $\Delta y_{i,t,t+1}=\Delta f(\theta_x,x_{i,t,t+1})+\Delta v_{i,t,t+1}-\Delta e_{i,t,t+1}$ in the spirit of \cite{chen2014consistent,belotti2012consistent}. Future research can extend this methods to the cases of persistent and transitory shocks from \cite{colombi2014closed} and \cite{kumbhakar2014technical}.
\section{Conclusions}
In this paper we study sensitivity analysis of Stochastic Frontier Models. We elaborate relaxations of the baseline assumptions in the Stochastic Frontier Models and characterize the identified set under this relaxations. Furthermore, we derive the breakdown frontier for a relevant parameter of interest, the average inefficiency of a production unit. We show an application of the procedures on a well known dataset, and make the code available for the interested practitioner.
\par
For a particular average level of inefficiency we find a downward slopping curve which exhibits the tradeoffs between the relaxations of the distributional assumption of the efficiency $e$ and of the error term $v$, that maintain the conclusions. The nonlinearity of the decreasing curve means that for a given relaxation of one assumption, say on $b$, it is required a reduction in the violation of the other ($c$) for the conclusions to hold, that varies depending on the value of $b$. Its convexity implies that a relaxation of $b$ when we allowed a small violation of that assumption permits a sizable reduction in the violation of the other assumption. And conversely, the same relaxation of $b$ when it is already broadly relaxed only allows for a small gain in the reduction of the violation of $c$.
\bibliographystyle{chicago}
\bibliography{references_all}