R/model_cs_tnorm_positive.R

In vh-d/SFAt: SFAt: Stochastic frontier analysis on cross-section and panel data

# NORMAL / T-NORMAL MODEL- CROSS SECTION DATA ----------------------------

par_cs_tnormp_check <- function(params,
                                indeces,
                                y, X,
                                CM = NULL,
                                CV_u,
                                CV_v) {

  # extract parameters from parameter vector
  n_f_coeff    <- ncol(X) # number of coeffs for frontier model
  n_cv_u_coeff <- if (is.matrix(CV_u)) ncol(CV_u) else length(CV_u)  # number of coeffs for conditional variance of inefficiency term model
  n_cv_v_coeff <- if (is.matrix(CV_v)) ncol(CV_v) else length(CV_v) # number of coeffs for conditional variance of the symmetric error model
  n_cm_coeff   <- if (is.matrix(CM)) ncol(CM) else length(CM) # number of coeffs for conditional mean model

  if (length(params) != n_f_coeff + n_cm_coeff + n_cv_u_coeff + n_cv_v_coeff) {
    stop("Incorrect nuber of parameters. ",
         n_f_coeff, "+", n_cm_coeff, "+", n_cv_u_coeff, "+", n_cv_v_coeff, " needed, but ", length(params), " supplied.")
  }

  return(TRUE)
}

t_par_cs_tnormp <- function(pars){
  # pars <- sqrt(exp(pars))
  # names(pars) <- c("sigma_v")
  # return(pars)
  return(NULL)
}


# likelihood --------------------------------------------------------------

# Advanced version of (Battese and Coelli, 1995) and (Huang and Liu, 1994) models
# heterogeneity in efficiency term: endogeneous location parameter mu
# implemented as Hadri et al. 2003
# parameters: f_coeff, cm_coeff, cv_u_coeff, cv_v_coeff

ll_cs_tnormp_contrib <- function(params,
                                 indeces,
                                 y, X,
                                 CV_u,
                                 CV_v,
                                 CM,
                                 ineff,
                                 deb = F) {

  if (deb) {
    cat("Parameters: ", params)
  }

  f_coeff    <- params[              1: indeces[1]]
  cv_u_coeff <- params[(indeces[1] + 1):indeces[2]]
  cv_v_coeff <- params[(indeces[2] + 1):indeces[3]]
  cm_coeff   <- params[(indeces[3] + 1):indeces[4]]

  sigma_u <- as.vector(exp(CV_u %*% cv_u_coeff))
  sigma_v <- as.vector(exp(CV_v %*% cv_v_coeff))

  sigma2_u <- sigma_u^2
  sigma2_v <- sigma_v^2

  sigma2 <- sigma2_v + sigma2_u # variance of composite error term
  sigma <- sqrt(sigma2)
  sigma_ast <- sigma_u * sigma_v / sigma

  if (deb) cat("Total of ", length(params), " parameters: \n",
               "frontier coeffs: ", paste(f_coeff), "\n",
               "cm coeffs: ",       paste(cm_coeff), "\n",
               "cv_u coeffs: ",     paste(cv_u_coeff), "\n",
               "cv_v coeffs: ",     paste(cv_v_coeff), "\n",
               "sigma_u: ", sigma_u, "\n",
               "sigma_v: ", sigma_v, "\n")

  N <- length(y)

  Zdelta <- as.vector(exp(CM %*% cm_coeff)) # fitted means of inefficiency term

  eps <- -ineff*as.vector(y - (X %*% f_coeff)) # composite error terms

  mu_ast <- (sigma2_v*Zdelta - sigma2_u*eps) / sigma2

  if (deb) {
    cat(length(Zdelta)    == N, "\n",
        length(eps)       == N, "\n",
        length(mu_ast)    == N, "\n",
        length(sigma_ast) == N, "\n",
        length(sigma2)    == N, "\n",
        length(sigma_u)   == N, "\n")
  }

  lli <- -0.5*log(2*pi) - log(sigma) - 0.5 * ((eps + Zdelta)^2)/sigma2 +
    log(pnorm(mu_ast / sigma_ast)) - log(pnorm(Zdelta / sigma_u))

  return(lli)
}


ll_cs_tnormp <- function(params,
                         indeces,
                         y, X,
                         CV_u,
                         CV_v,
                         CM,
                         ineff,
                         minmax = -1,
                         deb = F) {

  # compute loglik contributions
  lli <-
    ll_cs_tnormp_contrib(params,
                         indeces,
                         y, X,
                         CV_u,
                         CV_v,
                         CM,
                         ineff,
                         deb)

  if (deb) cat("Misbehaving loglikelihood contributions: ", "\n",
               which(!is.finite(lli)), "\n",
               lli[!is.finite(lli)],  "\n")

  if (deb) cat("Suspicious loglikelihood contributions: ", "\n",
               which(lli < quantile(lli, 0.05)), "\n",
               lli[lli < quantile(lli, 0.05)],  "\n")

  # sum to total loglikelihood
  ll <- sum(lli)
  if (deb) cat("Loglikelihood: ", ll,  "\n")

  return(minmax*ll) # return -loglikelihood for optimization of minimum
}



# gradient functions -----------------------------------------

# analyticaly derived gradient function
# todo
.g_cs_tnormp_analytic <- function(params,
                                  indeces,
                                  y, X,
                                  CV_u,
                                  CV_v,
                                  CM,
                                  ineff,
                                  minmax = -1,
                                  deb = F) {

  f_coeff    <- params[             1:indeces[1]]
  cv_u_coeff <- params[(indeces[1] + 1):indeces[2]]
  cv_v_coeff <- params[(indeces[2] + 1):indeces[3]]
  cm_coeff   <- params[(indeces[3] + 1):indeces[4]]

  sigma_u <- as.vector(exp(CV_u %*% cv_u_coeff))
  sigma_v <- as.vector(exp(CV_v %*% cv_v_coeff))

  sigma2_u <- sigma_u^2
  sigma2_v <- sigma_v^2

  sigma2 <- sigma2_v + sigma2_u # variance of composite error term
  sigma <- sqrt(sigma2)
  sigma_ast <- sigma_u * sigma_v / sigma

  N <- length(y)

  Zdelta <- as.vector(exp(CM %*% cm_coeff)) # fitted means of inefficiency term

  eps <- -ineff*as.vector(y - (X %*% f_coeff)) # composite error terms

  mu_ast <- (sigma2_v*Zdelta - sigma2_u*eps) / sigma2

  pdfnorm1 <- dnorm(mu_ast / sigma_ast)
  cdfnorm1 <- pnorm(mu_ast / sigma_ast)

  pdfnorm2 <- dnorm(Zdelta / sigma_u)
  cdfnorm2 <- pnorm(Zdelta / sigma_u)

  g_f_coeff <- -ineff*(pdfnorm1/cdfnorm1 * sigma_u/(sigma*sigma_v) + (eps + Zdelta)/sigma2) %*% X
  g_cv_u_coeff <- (-sigma2_u/sigma2 + (sigma2_u * (eps + Zdelta)^2) / (sigma2^2) + ((pdfnorm1/cdfnorm1) * (-2*sigma2_u*eps    - (sigma2_u/sigma2 + 1) * (sigma2_v*Zdelta - sigma2_u*eps))/(sigma*sigma_u*sigma_v)) + pdfnorm2/cdfnorm2 * Zdelta/sigma_u) %*% CV_u
  g_cv_v_coeff <- (-sigma2_v/sigma2 + (sigma2_v * (eps + Zdelta)^2) / (sigma2^2) + ((pdfnorm1/cdfnorm1) * (+2*sigma2_v*Zdelta - (sigma2_v/sigma2 + 1) * (sigma2_v*Zdelta - sigma2_u*eps))/(sigma*sigma_u*sigma_v))) %*% CV_v
  g_cm_coeff <-   (-1/sigma2 * (eps + Zdelta) + pdfnorm1/cdfnorm1 * sigma2_v / (sigma*sigma_u*sigma_v) - pdfnorm2/cdfnorm2 * 1/sigma_u) %*% CM

  # return(minmax*c(g_f_coeff, g_cv_u_coeff, g_cv_v_coeff, g_cm_coeff))
  return(minmax*c(g_f_coeff,
                  g_cv_u_coeff,
                  g_cv_v_coeff,
                  g_cm_coeff))
}

# finite-difference gradient function
g_cs_tnormp_fd <- function(params,
                           indeces,
                           y, X,
                           CV_u,
                           CV_v,
                           CM,
                           ineff,
                           minmax,
                           deb = F) {
  n <- length(params)
  hh <- matrix(0, n, n)
  diag(hh) <- .Machine$double.eps^(1/3)

  sapply(1:n, function(i) {
    (ll_cs_tnormp(params + hh[i, ], indeces, y, X, CV_u, CV_v, CM, ineff, minmax, deb) -
       ll_cs_tnormp(params - hh[i, ], indeces, y, X, CV_u, CV_v, CM, ineff, minmax, deb)) / (2 * .Machine$double.eps^(1/3))})
}

# until analytic is corrected
g_cs_tnormp_analytic <- g_cs_tnormp_fd


# efficiency --------------------------------------------------------------

u_cs_tnormp <- function(object, estimator) {
  # extract sigmas from the model object
  sigma_u <- as.vector(exp(object$data$CV_u %*% object$coeff_cv_u))
  sigma_v <- as.vector(exp(object$data$CV_v %*% object$coeff_cv_v))

  sigma2_u <- sigma_u^2
  sigma2_v <- sigma_v^2
  sigma2   <- sigma2_u + sigma2_v

  sigma_ast <- sqrt(sigma2_u * sigma2_v / sigma2)

  mu <- as.vector(exp(object$data$CM %*% object$coeff_cm))
  mu_ast <- (object$ineff * object$residuals * sigma2_u + sigma2_v * mu) / sigma2

  u <- switch(estimator,
              ME = pmax(mu_ast, 0),
              BC = -log(exp(-mu_ast + 0.5 * sigma_ast^2)*(pnorm(-sigma_ast + mu_ast/sigma_ast)/pnorm(mu_ast/sigma_ast))),
              JLMS = mu_ast + sigma_ast*(dnorm(mu_ast/sigma_ast)/pnorm(mu_ast/sigma_ast)),
              ... = stop(paste0("Unknown type ", estimator," of conditional mean estimator.")))

  return(u)
}