R/truncreg.R

In npsf: Nonparametric and Stochastic Efficiency and Productivity Analysis

truncreg <- function(formula, data, subset, 
                     ll = -Inf, ul = Inf, 
                     lmtol = .Machine$double.eps, maxiter = 150, 
                     marg.eff = FALSE, 
                     print.level = 1) {

 YX <- .prepareYXllul(formula = formula, ll = ll, ul = ul, data, subset, sysnframe = sys.nframe())
 y <- YX$Y
 x <- YX$X
 a <- YX$LL
 b <- YX$UL
 if(a[1] != -Inf & b[1] != Inf){
   model <- "ab"
 } else if (a[1] != -Inf) {
   model <- "a"
 } else if (b[1] != Inf) {
   model <- "b"
 }
 # sample size
 n <- YX$n
 n.full <- YX$n.full
 if(print.level >= 1){
  cat("\n(Note: ",n.full-n," (out of ",n.full,") observations truncated) \n\n", sep = "")
 }
 
 # initial values
 beta0 <- 0.999999*solve(t(x) %*% x) %*% t(x) %*% y
 sigma0 <- 1.000001*sqrt( sum((y - x %*% beta0)^2)/(n-ncol(x)) )
 names(sigma0) <- "/sigma"
 theta0 <- rbind( beta0, sigma0)
 K1 <- length(theta0)
 # initial LL
 XBeta     <- x %*% beta0
 yXBeta    <- y - XBeta
 if(model == "ab"){
   aXBeta  <- a - XBeta
   bXBeta  <- b - XBeta
   hdenom  <- pnorm(bXBeta/sigma0) - pnorm(aXBeta/sigma0)
 } else if (model == "a") {
   aXBeta  <- a - XBeta
   hdenom  <- 1 - pnorm(aXBeta/sigma0)
 } else if (model == "b") {
   bXBeta  <- b - XBeta
   hdenom  <- pnorm(bXBeta/sigma0)
 }
 l_init    <- sum( -0.5*log(2*pi*sigma0^2) - 0.5/sigma0^2*yXBeta^2 - log( hdenom ) )
#  # l_init <- .ll.trunc (y, x, beta0, sigma0, a, b)
 a0 <- a
 b0 <- b
 
 # print(summary(a))
 # print(summary(b))
 
 Delta <- 1
 mylm  <- 1
 step  <- -1
 time.05 <- proc.time()
 while ( abs(mylm) > lmtol & step <= maxiter){
  # print(step)
  step <- step + 1
  # print(step)
  # gradient
  if(model == "ab"){
    normbxb <- dnorm(bXBeta/sigma0)
    normaxb <- dnorm(aXBeta/sigma0)
    h0      <- (         normbxb -          normaxb) / hdenom
    h1      <- (  bXBeta*normbxb -   aXBeta*normaxb) / hdenom
    h2      <- (bXBeta^2*normbxb - aXBeta^2*normaxb) / hdenom
    h3      <- (bXBeta^3*normbxb - aXBeta^3*normaxb) / hdenom
  } else if (model == "a") {
    normaxb <- dnorm(aXBeta/sigma0)
    h0      <- (         -normaxb) / hdenom
    h1      <- (-aXBeta  *normaxb) / hdenom
    h2      <- (-aXBeta^2*normaxb) / hdenom
    h3      <- (-aXBeta^3*normaxb) / hdenom
  } else if (model == "b") {
    normbxb <- dnorm(bXBeta/sigma0)
    h0      <- (         normbxb) / hdenom
    h1      <- (  bXBeta*normbxb) / hdenom
    h2      <- (bXBeta^2*normbxb) / hdenom
    h3      <- (bXBeta^3*normbxb) / hdenom
  }
  # cat.print(h1)
  # h0 <- .h.trunc(al = aXBeta/sigma0, be = bXBeta/sigma0, a1 = 1,        b1 = 1,        a = a0, b = b0 )
  # h1 <- .h.trunc(al = aXBeta/sigma0, be = bXBeta/sigma0, a1 = aXBeta,   b1 = bXBeta,   a = a0, b = b0 )
  # h2 <- .h.trunc(al = aXBeta/sigma0, be = bXBeta/sigma0, a1 = aXBeta^2, b1 = bXBeta^2, a = a0, b = b0 )
  # h3 <- .h.trunc(al = aXBeta/sigma0, be = bXBeta/sigma0, a1 = aXBeta^3, b1 = bXBeta^3, a = a0, b = b0 )
  # h0 <- as.matrix(h0)
  # h1 <- as.matrix(h1)
  # h2 <- as.matrix(h2)
  # h3 <- as.matrix(h3)
  dl <- rbind(
   sigma0 ^ (-2) * t(x) %*% yXBeta + sigma0 ^ (-1) * t(x) %*% h0,
   -n * sigma0 ^ (-1) + sigma0 ^ (-3) * sum(yXBeta^2) + sigma0 ^
    (-2) * sum(h1)
  )
  # cat.print(dl)
  # KxN * Nx1
  # colnames(dl) <- "gradient"
  # print(dl)
  # Hessian
  # ddlA11 <- -sigma0^(-2)*t(x) %*% (diag( 1 - as.vector(sigma0^(-1) * h1 + h0^2) ) ) %*% x
  ddlA11 <- (-sigma0^(-2)*t(x * (1 - as.vector(sigma0^(-1) * h1 + h0^2) ) )  ) %*% x
  ddlA12 <- t(x) %*% (-2*sigma0^(-3)* yXBeta - sigma0^(-2)*h0 + sigma0^(-4)*h2 + sigma0^(-3)*h0*h1)
  # ddlA21 <- t(-2*sigma0^(-3)* yXBeta - sigma0^(-2)*h0 + sigma0^(-4)*h2 + sigma0^(-3)*h0*h1) %*% x
  ddlA22 <- n*sigma0^(-2) - 3*sigma0^(-4)* sum(yXBeta^2) - 2*sigma0^(-3) * sum(h1) + sigma0^(-5) * sum(h3) + sigma0^(-4) * sum(h1^2)
  ddl <- rbind( cbind(ddlA11, ddlA12), cbind(t(ddlA12), ddlA22) )
  # cat.print(ddl)
  iddl <- solve(ddl)
  theta0 <- theta0 - iddl %*% dl
  #print(theta0)
  beta0 <- theta0[-K1, , drop = FALSE]
  sigma0 <- theta0[K1]
  # sig0 <- sqrt(-diag(iddl))
  #print(cbind(dl, theta0))
  XBeta   <- x %*% beta0
  yXBeta  <- y - XBeta
  if(model == "ab"){
    aXBeta  <- a - XBeta
    bXBeta  <- b - XBeta
    hdenom  <- pnorm(bXBeta/sigma0) - pnorm(aXBeta/sigma0)
  } else if (model == "a") {
    aXBeta  <- a - XBeta
    hdenom  <- 1 - pnorm(aXBeta/sigma0)
  } else if (model == "b") {
    bXBeta  <- b - XBeta
    hdenom  <- pnorm(bXBeta/sigma0)
  }
  likelihood  <- sum( -0.5*log(2*pi*sigma0^2) - 0.5/sigma0^2*yXBeta^2 - log( hdenom ) )
  # likelihood <- .ll.trunc (y,x, beta0, sigma0,a,b)
  # print(likelihood, digits = 4)
  mylm <- as.vector( t(dl)%*%-iddl%*%dl )
  if ( is.na(mylm) ) stop("cannot compute an improvement")
  # cat(" Iteration ",step,": LM criterion = ",mylm,", LL = ",l_init," \n", sep = "")
  
  if ( print.level >= 1) {
   if ( print.level >= 2){
    cat(paste("Iteration ",formatC(step, width = 2)," (hessian is analytical, g*inv(H)*g' = ",formatC(mylm, format = "E", digits = 1),"): log likelihood = ",format(l_init, nsmall = 5),"\n", sep = ""), sep = "")
   } else {
    cat(paste("Iteration ",formatC(step, width = 2),": log likelihood = ",format(l_init, nsmall = 5),"\n", sep = ""), sep = "")
   }
  }
  Delta <- abs(likelihood - l_init)
  l_init <- likelihood
 } # the end of while
 
 time.06 <- proc.time()
 est.time.sec <- (time.06-time.05)[1]
 names(est.time.sec) <- "sec"
 if(print.level >= 1){
  .timing(est.time.sec, "\nLog likelihood maximization completed in\n")
  # cat("___________________________________________________\n")
 }
 
 if(step > maxiter) cat(" \n Not converged in ",maxiter," iterations \n\n", sep = "")
 # cat(" \n", sep = "")

 if ( print.level >= 1){
  cat("\nConvergence given g*inv(H)*g' = ",formatC(abs(mylm), format = "e", digits = 3)," < lmtol (",lmtol,")\n", sep = "")
  cat(paste("\nFinal log likelihood = ",formatC(l_init,digits=7,format="f"),"\n\n", sep = ""), sep = "")
  cat("Truncated regression for cross-sectional data\n", sep = "")
  cat("Limits:\n", sep = "")
  #  check if infinite
  ll.inf <- sum(is.infinite(YX$LL))
  ul.inf <- sum(is.infinite(YX$UL))
  
  if (ll.inf == n) {
   #  if all are infinite
   cat("lower limit for left-truncation  = ",YX$LL[1],"\n", sep = "")
  } else {
   if (sd(YX$LL[is.finite(YX$LL)]) != 0) {
    cat("lower limit for left-truncation:", sep = "")
    print(summary(YX$LL))
   } else {
    cat("lower limit for left-truncation  = ",YX$LL[1],"\n", sep = "")
   }
  }
  
  if (ul.inf == n) {
   #  if all are infinite
   cat("upper limit for right-truncation = ",YX$UL[1],"\n", sep = "")
  } else {
   if (sd(YX$UL[is.finite(YX$UL)]) != 0) {
    cat("upper limit for right-truncation:", sep = "")
    print(summary(YX$UL))
   } else {
    cat("upper limit for right-truncation = ",YX$UL[1],"\n", sep = "")
   }
  }
  
  cat("Number of observations (used in regression)               = ",n,"\n", sep = "")
  cat("Number of truncated observations (not used in regression) = ",n.full-n,"\n", sep = "")
  cat("\nEstimation results:\n\n", sep = "")
 }
 

 sig0 <- sqrt(-diag(iddl))

 # marginal effects
 # lmd <- NULL
 if(marg.eff){
  m.eff <- matrix(nrow = n, ncol = length(beta0))
  for(qq in 1:n)
  {
   x0 <- x[qq, ,drop = FALSE]
   a0 <- ifelse(length(a)==1, a, a[qq])
   b0 <- ifelse(length(b)==1, b, b[qq])

   if(b0 == Inf)
   {
    ai <- (a0 - x0 %*% beta0) / sigma0
    lambdai <- dnorm(ai) / (1 - pnorm(ai))
    # print(c(ai,lambdai))
    m.eff[qq,] <- as.matrix(1 - lambdai^2 + ai*lambdai) %*% t(beta0)
   } else if (a0 == -Inf)
   {
    ai <- (b0 - x0 %*% beta0) / sigma0
    lambdai <- -dnorm(ai) / pnorm(ai)
    m.eff[qq,] <- as.matrix(1 - lambdai^2 + ai*lambdai) %*% t(beta0)
   } else
   {
    ai1 <- (a0 - x0 %*% beta0) / sigma0
    ai2 <- (b0 - x0 %*% beta0) / sigma0
    lambdai <- (dnorm(ai2)-dnorm(ai1)) / (pnorm(ai2) - pnorm(ai1))
    m.eff[qq,] <- as.matrix(1 - lambdai^2 + ai2*lambdai - (b0-a0)*dnorm(ai1) / sigma0 / ( pnorm(ai2) - pnorm(ai1) )) %*% t(beta0)
   }
   # lmd <- c(lmd, lambdai)
  } # end of qq loop
  colnames(m.eff) <- rownames(beta0)
 }

 coeffs <- cbind(theta0, sig0,theta0/sig0,(1-pnorm(abs(theta0/sig0),lower.tail =TRUE))*2)
 rownames(coeffs)[length(beta0)+1] <- "/sigma"
 colnames(coeffs) <- c("Estimate","Std. Error","z","Pr(>|z|)")

 if (print.level >= 1){
  printCoefmat(coeffs)
 }

 tymch <- list(call = match.call(), model = "truncreg.cs", coef = coeffs[,1], table = coeffs, vcov = -iddl, ll = l_init, lmtol = lmtol, LM = mylm, esttime = est.time.sec, sigma = theta0[length(theta0)], LL = a, UL = b, n = n, n.full = n.full, nontruncsample = as.vector(YX$nontruncsample), esample = YX$esample)
 
 if(marg.eff){
  tymch$marg.effects <- m.eff
 }

 class(tymch) <- "npsf"
 return(tymch)
}