robustSFA (version 0.2.0)

rsfa: Robust Estimation of Stochastic Frontier Models

Description

This function performs the robust estimation for stochastic frontier models using the Minimum Density Power Divergence Estimator (MDPDE). The MDPDE is particularly useful when outliers in a dataset distort estimation results, especially regarding technical efficiencies. The parameter \(\alpha\) in the MDPDE plays a crucial role in mitigating the impact of outliers when estimating coefficients and technical efficiencies. It actually controls the trade-off between robustness and efficiency of the MDPDE. The current estimation process is based on the following model: $$Y=\beta_0 +\beta_1 X_1 +\cdots +\beta_p X_p + V - U:=g(X,\beta)+V-U,$$ where \(V\) follows the normal distribution with mean 0 and variance \(\sigma_v^2\), and \(U\) follows the half-normal distribution with variance \(\sigma_u^2\). The MDPDE with \(\alpha\) for the model above is given as follows: $$ \underset{\theta \in \Theta}{\operatorname{argmin}} \, \frac{1}{\sigma^\alpha} \bigg\{ \frac{\sqrt{2}}{\sqrt{\pi}} \int e^{-(1+\alpha)\frac{z^2}{2}}\Phi^{1+\alpha}(-z\lambda)dz$$ $$ \hspace{5cm} - \Big(1+\frac{1}{\alpha}\Big)\frac{1}{n} \sum_{i=1}^n e^{-\alpha\frac{(Y_i-g(X_i,\beta))^2}{2\sigma^2}} \Phi^\alpha\Big(-\frac{Y_i-g(X_i,\beta)}{\sigma}\lambda\Big) \bigg\}, $$ where \(\theta=(\beta_0, \beta_1, \cdots, \beta_p, \sigma^2, \lambda)\) with \(\sigma^2=\sigma^2_v+\sigma^2_u\) and \(\lambda=\sigma_u/\sigma_v\). For more details, refer to Song et al. (2017).

Usage

rsfa(formula, data = NULL, alpha = 0, se = TRUE)

Value

A list containing robust estimation results: estimates, standard errors, and residuals, as well as technical efficiencies.

call

The matched call.

coefficients

A named vector of parameter estimates.

se

A vector of estimated standard errors.

residuals

A vector of residuals.

efficiencies

A vector of estimated technical efficiencies calculated using the estimator of Battese and Coelli (1988).

sigma2_V

Numeric. Estimated value of \(\sigma^2_V\).

sigma2_U

Numeric. Estimated value of \(\sigma^2_U\).

alpha

Numeric. \(\alpha\) used in the MDPDE.

Arguments

formula

A symbolic description of the model to be estimated, specified using the standard R formula syntax (e.g., y ~ x1 + x2).

data

A data frame containing the variables in the model.

alpha

A numeric value controlling the robustness of the MDPDE. When alpha is zero (the default), 'rsfa' returns the ML estimates. It is not necessary for \(\alpha\) to be large. Values of \(\alpha\) close to zero are known to be efficient and sufficiently robust. The usually recommended range of \(\alpha\) is 0.1 to 0.3. A function for selecting the optimal \(\alpha\) is also provided in this package. For this purpose, see the optimal_alpha() function.

se

Logical. If TRUE (the default), standard errors are presented.

References

Battese, G.E. and Coelli, T. (1988). Prediction of firm-level technical efficiencies with a generalized frontier production function and panel data. Journal of Econometrics, 38, 387-399.

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

Song, J., Oh, D-h., and Kang, J. (2017). Robust estimation in stochastic frontier models. Computational Statistics & Data Analysis, 105, 243-267.

Examples

# \donttest{
## Example using the 'riceProdPhil' dataset from the `frontier` package
library(frontier)
data(riceProdPhil)
summary(riceProdPhil)


## ML estimates and MDPD estimates with alpha=0.1
my.model <- log(PROD) ~ log(AREA) + log(LABOR) + log(NPK) + log(OTHER)
fit.ml <- rsfa(my.model, data = riceProdPhil)
summary(fit.ml)

fit.mdpde<- rsfa(my.model, data = riceProdPhil, alpha = 0.1)
summary(fit.mdpde)


## Comparison of the technical efficiencies from MLE and MDPDE
eff.ml<-efficiency(fit.ml)
eff.mdpde<-efficiency(fit.mdpde)
plot(eff.ml, eff.mdpde, pch=20,
     xlab="efficiency from MLE", ylab="efficiency from MDPDE")
abline(0,1, lty=2)


## Data with a single outlying observation
riceProdPhil2 <- riceProdPhil
riceProdPhil3 <- riceProdPhil
idx <- which.max(riceProdPhil$PROD)
riceProdPhil2$PROD[idx] <- riceProdPhil$PROD[idx]*10
riceProdPhil3$PROD[idx] <- riceProdPhil$PROD[idx]/100


## The technical efficiencies for "riceProdPhil2" data
fit.ml2<- rsfa(my.model, data = riceProdPhil2)
eff.ml2 <- efficiency(fit.ml2)
fit.mdpde2<- rsfa(my.model, data = riceProdPhil2, alpha = 0.1)
eff.mdpde2 <- efficiency(fit.mdpde2)

plot(eff.ml, eff.ml2, xlim = c(0, 1), ylim = c(0, 1), col="red", pch=20,
     xlab = "efficiency from MLE (riceProdPhil)",
     ylab = "efficiency (riceProdPhil2)")
points(eff.ml, eff.mdpde2, col="blue", pch=20)
abline(0,1,lty=2)
legend("topleft", legend = c("MLE", "MDPDE"), col = c("red", "blue"), pch = 20, bty="n")     
title(main =
          "Efficiency comparison: MLE(no outlier) vs MLE and MPDPE (one over-performing producer)")


## The technical efficiencies for "riceProdPhil3" data
fit.ml3<- rsfa(my.model, data = riceProdPhil3)
eff.ml3 <- efficiency(fit.ml3)
fit.mdpde3<- rsfa(my.model, data = riceProdPhil3, alpha = 0.1)
eff.mdpde3 <- efficiency(fit.mdpde3)

plot(eff.ml, eff.ml3, xlim = c(0, 1), ylim = c(0, 1), col="red", pch=20,
     xlab = "efficiency from MLE (riceProdPhil)",
     ylab = "efficiency (riceProdPhil3)")
points(eff.ml, eff.mdpde3, col="blue", pch=20)
abline(0,1,lty=2)
legend("topleft", legend = c("MLE", "MDPDE"), col = c("red", "blue"), pch = 20, bty="n")     
title(main =
          "Efficiency comparison: MLE(no outlier) vs MLE and MPDPE (one under-performing producer)")
# }