This method fits non-parametric stochastic frontier models. The data-generating process
is assumed to be of the form
$$\ln y_i = \ln f(x_i) + v_i - u_i,$$
where \(y_i\) is the \(i\)th observation of output, \(f\) is a continuous
function, \(x_i\) is the \(i\)th observation of input, \(v_i\) is a
normally-distributed error term (\(v_i\sim N(0, \sigma_v^2)\)), and \(u_i\) is a
normally-distributed error term truncated below at zero (\(u_i\sim N^+(0, \sigma_u)\)).
Aigner et al. developed methods to decompose
\(\varepsilon_i = v_i - u_i\) into its basic components.
This procedure first fits the mean of the data using fit.mean,
producing estimates of output \(\hat{y}\). Log-proportional errors are calculated as
$$\varepsilon_i = \ln(y_i / \hat{y}_i).$$
Following Aigner et al. (1977), parameters of one- and two-sided error distributions
are estimated via maximum likelihood. First,
$$\hat{\sigma}^2 = \frac1N \sum_{i=1}^N \varepsilon_i^2.$$
Then, \(\hat{\lambda}\) is estimated by solving
$$\frac1{\hat{\sigma}^2} \sum_{i=1}^N \varepsilon_i\hat{y}_i + \frac{\hat{\lambda}}{\hat{\sigma}} \sum_{i=1}^N \frac{f_i^*}{1 - F_i^*}y_i = 0,$$
where \(f_i^*\) and \(F_i^*\) are standard normal density and distribution function,
respectively, evaluated at \(\varepsilon_i\hat{\lambda}\hat{\sigma}^{-1}\). Parameters of
the one- and two-sided distributions are found by solving the identities
$$\sigma^2 = \sigma_u^2 + \sigma_v^2$$
$$\lambda = \frac{\sigma_u}{\sigma_v}.$$
Mean efficiency over the sample is given by
$$\exp\left(-\frac{\sqrt{2}}{\sqrt{\pi}}\right)\sigma_u,$$
and modal efficiency for each observation is given by
$$-\varepsilon(\sigma_u^2/\sigma^2).$$