R/efficiencies.R

In sfaR: Stochastic Frontier Analysis Routines

################################################################################
#                                                                              #
# R functions for the sfaR package                                             #
#                                                                              #
################################################################################

#------------------------------------------------------------------------------#
# Efficiency/Inefficiency estimation                                           #
# Models: + Cross sectional & Pooled data                                      #
#           -Stochastic Frontier Analysis                                      #
#           -Latent Class Stochastic Frontier Analysis                         #
#           -Sample selection correction for Stochastic Frontier Model         #
# Data: Cross sectional data & Pooled data                                     #
#------------------------------------------------------------------------------#

#' Compute conditional (in-)efficiency estimates of stochastic frontier models
#' 
#' \code{\link{efficiencies}} returns (in-)efficiency estimates of models 
#' estimated with \code{\link{sfacross}}, \code{\link{sfalcmcross}}, or 
#' \code{\link{sfaselectioncross}}.
#' 
#' @name efficiencies
#' 
#' @details In general, the conditional inefficiency is obtained following 
#' Jondrow \emph{et al.} (1982) and the conditional efficiency is computed 
#' following Battese and Coelli (1988). In some cases the conditional mode is 
#' also returned (Jondrow \emph{et al.} 1982). The confidence interval is 
#' computed following Horrace and Schmidt (1996), Hjalmarsson \emph{et al.} 
#' (1996), or Berra and Sharma (1999) (see \sQuote{Value} section).
#'
#' In the case of the half normal distribution for the one-sided error term,
#' the formulae are as follows (for notations, see the \sQuote{Details} section
#' of \code{\link{sfacross}} or \code{\link{sfalcmcross}}):
#'
#' \itemize{ \item The conditional inefficiency is: }
#' 
#' \deqn{E\left\lbrack u_i|\epsilon_i\right
#' \rbrack=\mu_{i\ast} + \sigma_\ast\frac{\phi
#' \left(\frac{\mu_{i\ast}}{\sigma_\ast}\right)}{
#' \Phi\left(\frac{\mu_{i\ast}}{\sigma_\ast}\right)}}
#'
#' where
#' 
#' \deqn{\mu_{i\ast}=\frac{-S\epsilon_i\sigma_u^2}{ \sigma_u^2 + \sigma_v^2}}
#'
#' and
#' 
#' \deqn{\sigma_\ast^2 = \frac{\sigma_u^2 \sigma_v^2}{\sigma_u^2 + \sigma_v^2}}
#' 
#' \itemize{ \item The Battese and Coelli (1988) conditional efficiency is
#' obtained with: } 
#' 
#' \deqn{E\left\lbrack\exp{\left(-u_i\right)}
#' |\epsilon_i\right\rbrack = \exp{\left(-\mu_{i\ast}+
#' \frac{1}{2}\sigma_\ast^2\right)}\frac{\Phi\left(
#' \frac{\mu_{i\ast}}{\sigma_\ast}-\sigma_\ast\right)}{
#' \Phi\left(\frac{\mu_{i\ast}}{\sigma_\ast}\right)}}
#' 
#' \itemize{ \item The reciprocal of the Battese and Coelli (1988) conditional 
#' efficiency is obtained with: } 
#' 
#' \deqn{E\left\lbrack\exp{\left(u_i\right)}
#' |\epsilon_i\right\rbrack = \exp{\left(\mu_{i\ast}+
#' \frac{1}{2}\sigma_\ast^2\right)} \frac{\Phi\left(
#' \frac{\mu_{i\ast}}{\sigma_\ast}+\sigma_\ast\right)}{
#' \Phi\left(\frac{\mu_{i\ast}}{\sigma_\ast}\right)}}
#' 
#' \itemize{ \item The conditional mode is computed using: }
#' 
#' \deqn{M\left\lbrack u_i|\epsilon_i\right
#' \rbrack= \mu_{i\ast} \quad \hbox{For} \quad 
#' \mu_{i\ast} > 0}
#' 
#' and
#' 
#' \deqn{M\left\lbrack u_i|\epsilon_i\right
#' \rbrack= 0 \quad \hbox{For} \quad \mu_{i\ast} \leq 0}
#' 
#' \itemize{ \item The confidence intervals are obtained with: }
#' 
#' \deqn{\mu_{i\ast} + I_L\sigma_\ast \leq 
#' E\left\lbrack u_i|\epsilon_i\right\rbrack \leq 
#' \mu_{i\ast} + I_U\sigma_\ast }
#'
#' with \eqn{LB_i = \mu_{i*} + I_L\sigma_*} and 
#' \eqn{UB_i = \mu_{i*} + I_U\sigma_*}
#'
#' and
#' 
#' \deqn{I_L = \Phi^{-1}\left\lbrace 1 -
#' \left(1-\frac{\alpha}{2}\right)\left\lbrack 1-
#' \Phi\left(-\frac{\mu_{i\ast}}{\sigma_\ast}\right)
#' \right\rbrack\right\rbrace }
#'
#' and
#' 
#' \deqn{I_U = \Phi^{-1}\left\lbrace 1-
#' \frac{\alpha}{2}\left\lbrack 1-\Phi
#' \left(-\frac{\mu_{i\ast}}{\sigma_\ast}\right)
#' \right\rbrack\right\rbrace}
#'
#' Thus
#' 
#' \deqn{\exp{\left(-UB_i\right)} \leq E\left
#' \lbrack\exp{\left(-u_i\right)}|\epsilon_i\right\rbrack 
#' \leq\exp{\left(-LB_i\right)}}
#' 
#' In the case of the sample selection, as underlined in Greene (2010), the 
#' conditional inefficiency could be computed using Jondrow \emph{et al.} (1982).
#' However, here the conditionanl (in)efficiency is obtained using the properties
#' of the closed skew-normal (CSN) distribution (Lai, 2015). The conditional
#' efficiency can be obtained using the moment generating functions of a CSN 
#' distribution (see Gonzalez-Farias \emph{et al.} (2004)). We have:
#' 
#' \deqn{E\left\lbrack\exp{\left(tu_i\right)}
#' |\epsilon_i\right\rbrack = M_{u|\epsilon}(t)=\frac{\Phi_2\left(\tilde{\mathbf{D}}
#' \tilde{\bm{\Sigma}}t; \tilde{\bm{\kappa}}, \tilde{\bm{\Delta}} + 
#' \tilde{\mathbf{D}}\tilde{\bm{\Sigma}}\tilde{\mathbf{D}}' \right)}{
#' \Phi_2\left(\mathbf{0}; \tilde{\bm{\kappa}}, \tilde{\bm{\Delta}} + 
#' \tilde{\mathbf{D}}\tilde{\bm{\Sigma}}\tilde{\mathbf{D}}'\right)}\exp{
#' \left(t\tilde{\bm{\pi}} + \frac{1}{2}t^2\tilde{\bm{\Sigma}}\right)}}
#' 
#' where \eqn{\tilde{\bm{\pi}} = \frac{-S\epsilon_i\sigma_u^2}{\sigma_v^2 + \sigma_u^2}}, 
#' \eqn{\tilde{\bm{\Sigma}} = \frac{\sigma_v^2\sigma_u^2}{\sigma_v^2 + \sigma_u^2}}, 
#' \eqn{\tilde{\mathbf{D}} = \begin{pmatrix} \frac{S\rho}{\sigma_v} \\ 1 \end{pmatrix}}, 
#' \eqn{\tilde{\bm{\kappa}} = \begin{pmatrix} - \mathbf{Z}'_{si}\bm{\gamma} - 
#' \frac{\rho\sigma_v\epsilon_i}{\sigma_v^2 + \sigma_u^2}\\ 
#' \frac{S\sigma_u^2\epsilon_i}{\sigma_v^2 + \sigma_u^2} \end{pmatrix}}, 
#' \eqn{\tilde{\bm{\Delta}} = \begin{pmatrix}1-\rho^2 & 0 \\ 0 & 0\end{pmatrix}}.
#' 
#' The derivation of the efficiency and the reciprocal efficiency is obtained by replacing
#' \eqn{t = -1} and \eqn{t =1}, respectively. To obtain the inefficiency as 
#' \eqn{E\left[u_i|\epsilon_i\right]} is more complicated as it requires the 
#' derivation of a multivariate normal cdf. We have:
#' 
#' \deqn{E\left[u_i|\epsilon_i\right] = \left. \frac{\partial M_{u|\epsilon}(t)}{\partial t}\right\rvert_{t = 0}}
#' 
#' Then
#' 
#' \deqn{E\left[u_i|\epsilon_i\right] = \tilde{\bm{\pi}} + 
#' \left(\tilde{\mathbf{D}}\tilde{\bm{\Sigma}}\right)'\frac{\Phi_2^*
#' \left(\mathbf{0}; \tilde{\bm{\kappa}}, \ddot{\bm{\Delta}}\right)}{
#' \Phi_2\left(\mathbf{0}; \tilde{\bm{\kappa}}, \ddot{\bm{\Delta}}\right)}}
#' 
#' where \eqn{\Phi_2^* \left(\mathbf{s}; \tilde{\bm{\kappa}}, \ddot{\bm{\Delta}}\right)=
#' \frac{\partial \Phi_2\left(\mathbf{s}; \tilde{\bm{\kappa}}, \ddot{\bm{\Delta}} \right)}{\partial \mathbf{s}}}
#'
#' @param object A stochastic frontier model returned
#' by \code{\link{sfacross}}, \code{\link{sfalcmcross}}, or 
#' \code{\link{sfaselectioncross}}.
#' @param level A number between between 0 and 0.9999 used for the computation
#' of (in-)efficiency confidence intervals (defaut = \code{0.95}). Only used
#' when \code{udist} = \code{'hnormal'}, \code{'exponential'}, \code{'tnormal'}
#' or \code{'uniform'} in \code{\link{sfacross}}.
#' @param newData Optional data frame that is used to calculate the efficiency 
#' estimates. If NULL (the default), the efficiency estimates are calculated 
#' for the observations that were used in the estimation. In the case of object of 
#' class \code{sfaselectioncross} 
#' @param ... Currently ignored.
#'
#' @return A data frame that contains individual (in-)efficiency estimates.
#' These are ordered in the same way as the corresponding observations in the
#' dataset used for the estimation.
#' 
#' \bold{- For object of class \code{'sfacross'} the following elements are returned:}
#'
#' \item{u}{Conditional inefficiency. In the case argument \code{udist} of
#' \link{sfacross} is set to \code{'uniform'}, two conditional inefficiency
#' estimates are returned: \code{u1} for the classic conditional inefficiency
#' following Jondrow \emph{et al.} (1982), and \code{u2} which is obtained when
#' \eqn{\theta/\sigma_v \longrightarrow \infty} (see Nguyen, 2010).}
#'
#' \item{uLB}{Lower bound for conditional inefficiency. Only when the argument
#' \code{udist} of \link{sfacross} is set to \code{'hnormal'},
#' \code{'exponential'}, \code{'tnormal'} or \code{'uniform'}.}
#'
#' \item{uUB}{Upper bound for conditional inefficiency. Only when the argument
#' \code{udist} of \link{sfacross} is set to \code{'hnormal'},
#' \code{'exponential'}, \code{'tnormal'} or \code{'uniform'}.}
#'
#' \item{teJLMS}{\eqn{\exp{(-E[u|\epsilon])}}. When the argument \code{udist} of
#' \link{sfacross} is set to \code{'uniform'}, \code{teJLMS1} =
#' \eqn{\exp{(-E[u_1|\epsilon])}} and \code{teJLMS2} = 
#' \eqn{\exp{(-E[u_2|\epsilon])}}. Only when \code{logDepVar = TRUE}.}
#'
#' \item{m}{Conditional model. Only when the argument \code{udist} of
#' \link{sfacross} is set to \code{'hnormal'}, \code{'exponential'},
#' \code{'tnormal'}, or \code{'rayleigh'}.}
#'
#' \item{teMO}{\eqn{\exp{(-m)}}. Only when, in the function \link{sfacross},
#' \code{logDepVar = TRUE} and \code{udist = 'hnormal'}, \code{'exponential'},
#' \code{'tnormal'}, \code{'uniform'}, or \code{'rayleigh'}.}
#'
#' \item{teBC}{Battese and Coelli (1988) conditional efficiency. Only when, in
#' the function \link{sfacross},  \code{logDepVar = TRUE}. 
#' In the case \code{udist = 'uniform'}, two 
#' conditional efficiency estimates are returned:
#' \code{teBC1} which is the classic conditional efficiency following 
#' Battese and Coelli (1988) and \code{teBC2} when 
#' \eqn{\theta/\sigma_v \longrightarrow \infty} (see Nguyen, 2010).}
#' 
#' \item{teBC_reciprocal}{Reciprocal of Battese and Coelli (1988) conditional 
#' efficiency. Similar to \code{teBC} except that it is computed as 
#' \eqn{E\left[\exp{(u)}|\epsilon\right]}.}
#'
#' \item{teBCLB}{Lower bound for Battese and Coelli (1988) conditional
#' efficiency. Only when, in the function \link{sfacross}, \code{logDepVar = TRUE} and 
#' \code{udist = 'hnormal'}, \code{'exponential'}, \code{'tnormal'},
#' or \code{'uniform'}.}
#'
#' \item{teBCUB}{Upper bound for Battese and Coelli (1988) conditional
#' efficiency. Only when, in the function \link{sfacross}, \code{logDepVar = TRUE} and 
#' \code{udist = 'hnormal'}, \code{'exponential'}, \code{'tnormal'},
#' or \code{'uniform'}.}
#' 
#' \item{theta}{In the case \code{udist = 'uniform'}. \eqn{u \in [0, \theta]}.}
#' 
#' \bold{- For object of class \code{'sfalcmcross'} the following elements are returned:}
#'
#' \item{Group_c}{Most probable class for each observation.}
#'
#' \item{PosteriorProb_c}{Highest posterior probability.}
#' 
#' \item{u_c}{Conditional inefficiency of the most probable class given the
#' posterior probability.}
#' 
#' \item{teJLMS_c}{\eqn{\exp{(-E[u_c|\epsilon_c])}}. Only when, in the function
#' \link{sfalcmcross} \code{logDepVar = TRUE}.}
#' 
#' \item{teBC_c}{\eqn{E\left[\exp{(-u_c)}|\epsilon_c\right]}. Only when, in the 
#' function \link{sfalcmcross} \code{logDepVar = TRUE}.}
#' 
#' \item{teBC_reciprocal_c}{\eqn{E\left[\exp{(u_c)}|\epsilon_c\right]}. Only 
#' when, in the function \link{sfalcmcross} \code{logDepVar = TRUE}.}
#'
#' \item{PosteriorProb_c#}{Posterior probability of class #.}
#'
#' \item{PriorProb_c#}{Prior probability of class #.}
#'
#' \item{u_c#}{Conditional inefficiency associated to class #, regardless of
#' \code{Group_c}.}
#' 
#' \item{teBC_c#}{Conditional efficiency 
#' (\eqn{E\left[\exp{(-u_c)}|\epsilon_c\right]}) associated to class #, 
#' regardless of \code{Group_c}. Only when, in the function
#' \link{sfalcmcross} \code{logDepVar = TRUE}.}
#' 
#' \item{teBC_reciprocal_c#}{Reciprocal conditional efficiency 
#' (\eqn{E\left[\exp{(u_c)}|\epsilon_c\right]}) associated to class #, 
#' regardless of \code{Group_c}. Only when, in the function
#' \link{sfalcmcross} \code{logDepVar = TRUE}.}
#'
#' \item{ineff_c#}{Conditional inefficiency (\code{u_c}) for observations in
#' class # only.}
#' 
#' \item{effBC_c#}{Conditional efficiency (\code{teBC_c}) for observations in
#' class # only.}
#' 
#' \item{ReffBC_c#}{Reciprocal conditional efficiency (\code{teBC_reciprocal_c}) 
#' for observations in class # only.}
#' 
#' \item{theta_c#}{In the case \code{udist = 'uniform'}. \eqn{u \in [0, \theta_{c\#}]}.}
#' 
#' \bold{- For object of class \code{'sfaselectioncross'} the following elements are returned:}
#' 
#' \item{u}{Conditional inefficiency.}
#'
#' \item{teJLMS}{\eqn{\exp{(-E[u|\epsilon])}}. Only when \code{logDepVar = TRUE}.}
#'
#' \item{teBC}{Battese and Coelli (1988) conditional efficiency. Only when, in
#' the function \link{sfaselectioncross}, 
#' \code{logDepVar = TRUE}.}
#' 
#' \item{teBC_reciprocal}{Reciprocal of Battese and Coelli (1988) conditional 
#' efficiency. Similar to \code{teBC} except that it is computed as 
#' \eqn{E\left[\exp{(u)}|\epsilon\right]}.}
#' 
# @author K Hervé Dakpo
#'
#' @seealso \code{\link{sfalcmcross}}, for the latent class stochastic frontier analysis
#' model fitting function using cross-sectional or pooled data.
#' 
#' \code{\link{sfacross}}, for the stochastic frontier analysis model
#' fitting function using cross-sectional or pooled data.
#' 
#' \code{\link{sfaselectioncross}} for sample selection in stochastic frontier 
#' model fitting function using cross-sectional or pooled data.
#'
#' @references Battese, G.E., and T.J. Coelli. 1988. Prediction of firm-level
#' technical efficiencies with a generalized frontier production function and
#' panel data. \emph{Journal of Econometrics}, \bold{38}:387--399.
#'
#' Bera, A.K., and S.C. Sharma. 1999. Estimating production uncertainty in
#' stochastic frontier production function models. \emph{Journal of
#' Productivity Analysis}, \bold{12}:187-210.
#' 
#' Gonzalez-Farias, G., Dominguez-Molina, A., Gupta, A. K., 2004. Additive 
#' properties of skew normal random vectors. 
#' \emph{Journal of Statistical Planning and Inference}. \bold{126}: 521-534.
#' 
#' Greene, W., 2010. A stochastic frontier model with correction 
#' for sample selection. \emph{Journal of Productivity Analysis}. \bold{34}, 15--24.
#'
#' Hjalmarsson, L., S.C. Kumbhakar, and A. Heshmati. 1996. DEA, DFA and SFA: A
#' comparison. \emph{Journal of Productivity Analysis}, \bold{7}:303-327.
#'
#' Horrace, W.C., and P. Schmidt. 1996. Confidence statements for efficiency
#' estimates from stochastic frontier models. \emph{Journal of Productivity
#' Analysis}, \bold{7}:257-282.
#'
#' Jondrow, J., C.A.K. Lovell, I.S. Materov, and P. Schmidt. 1982. On the
#' estimation of technical inefficiency in the stochastic frontier production
#' function model. \emph{Journal of Econometrics}, \bold{19}:233--238.
#' 
#' Lai, H. P., 2015. Maximum likelihood estimation of the stochastic frontier 
#' model with endogenous switching or sample selection. 
#' \emph{Journal of Productivity Analysis}, \bold{43}: 105-117.
#'
#' Nguyen, N.B. 2010. Estimation of technical efficiency in stochastic frontier
#' analysis. PhD Dissertation, Bowling Green State University, August.
#'
#' @examples
#' 
#' \dontrun{
#' ## Using data on fossil fuel fired steam electric power generation plants in the U.S.
#' # Translog SFA (cost function) truncated normal with scaling property
#' tl_u_ts <- sfacross(formula = log(tc/wf) ~ log(y) + I(1/2 * (log(y))^2) + log(wl/wf) +
#' log(wk/wf) + I(1/2 * (log(wl/wf))^2) + I(1/2 * (log(wk/wf))^2) + I(log(wl/wf) *
#' log(wk/wf)) + I(log(y) * log(wl/wf)) + I(log(y) * log(wk/wf)), udist = 'tnormal',
#' muhet = ~ regu, uhet = ~ regu, data = utility, S = -1, scaling = TRUE, method = 'mla')
#' eff.tl_u_ts <- efficiencies(tl_u_ts)
#' head(eff.tl_u_ts)
#' summary(eff.tl_u_ts)
#' }
#' 
#' @aliases efficiencies.sfacross
#' @export
#' @export efficiencies
# conditional efficiencies sfacross ----------
efficiencies.sfacross <- function(object, level = 0.95, newData = NULL,
  ...) {
  if (level < 0 || level > 0.9999) {
    stop("'level' must be between 0 and 0.9999", call. = FALSE)
  }
  if (!is.null(newData)) {
    if (!is.data.frame(newData)) {
      stop("argument 'newData' must be of class data.frame")
    }
    object$dataTable <- newData
    object$Nobs <- dim(newData)[1]
  }
  if (object$udist == "hnormal") {
    EffRes <- chalfnormeff(object = object, level = level)
  } else {
    if (object$udist == "exponential") {
      EffRes <- cexponormeff(object = object, level = level)
    } else {
      if (object$udist == "gamma") {
        EffRes <- cgammanormeff(object = object, level = level)
      } else {
        if (object$udist == "rayleigh") {
          EffRes <- craynormeff(object = object, level = level)
        } else {
          if (object$udist == "uniform") {
          EffRes <- cuninormeff(object = object, level = level)
          } else {
          if (object$udist == "tnormal") {
            if (object$scaling) {
            EffRes <- ctruncnormscaleff(object = object,
              level = level)
            } else {
            EffRes <- ctruncnormeff(object = object,
              level = level)
            }
          } else {
            if (object$udist == "lognormal") {
            EffRes <- clognormeff(object = object,
              level = level)
            } else {
            if (object$udist == "genexponential") {
              EffRes <- cgenexponormeff(object = object,
              level = level)
            } else {
              if (object$udist == "tslaplace") {
              EffRes <- ctslnormeff(object = object,
                level = level)
              } else {
              if (object$udist == "weibull") {
                EffRes <- cweibullnormeff(object = object,
                level = level)
              }
              }
            }
            }
          }
          }
        }
      }
    }
  }
  return(EffRes)
}

# conditional efficiencies sfalcmcross ----------
#' @rdname efficiencies
#' @aliases efficiencies.sfalcmcross
#' @export
efficiencies.sfalcmcross <- function(object, level = 0.95, newData = NULL,
  ...) {
  if (level < 0 || level > 0.9999) {
    stop("'level' must be between 0 and 0.9999", call. = FALSE)
  }
  if (!is.null(newData)) {
    if (!is.data.frame(newData)) {
      stop("argument 'newData' must be of class data.frame")
    }
    object$dataTable <- newData
    object$Nobs <- dim(newData)[1]
  }
  if (object$nClasses == 2) {
    EffRes <- cLCM2Chalfnormeff(object = object, level = level)
  } else {
    if (object$nClasses == 3) {
      EffRes <- cLCM3Chalfnormeff(object = object, level = level)
    } else {
      if (object$nClasses == 4) {
        EffRes <- cLCM4Chalfnormeff(object = object,
          level = level)
      } else {
        if (object$nClasses == 5) {
          EffRes <- cLCM5Chalfnormeff(object = object,
          level = level)
        }
      }
    }
  }
  return(EffRes)
}

# conditional efficiencies sfaselectioncross ----------
#' @rdname efficiencies
#' @aliases efficiencies.sfaselectioncross
#' @export
efficiencies.sfaselectioncross <- function(object, level = 0.95,
  newData = NULL, ...) {
  if (level < 0 || level > 0.9999) {
    stop("'level' must be between 0 and 0.9999", call. = FALSE)
  }
  if (!is.null(newData)) {
    if (!is.data.frame(newData)) {
      stop("argument 'newData' must be of class data.frame")
    }
    object$dataTable <- newData
    object$Ninit <- dim(newData)[1]
  }
  EffRes <- chalfnormeff_ss(object = object, level = level)
  return(EffRes)
}