rsfa {robustSFA} | R Documentation |
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)
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 |
se |
Logical. If TRUE (the default), standard errors are presented. |
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 |
sigma2_U |
Numeric. Estimated value of |
alpha |
Numeric. |
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
## 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)")