R/iqr_milp.R

In ChenyueLiu/ivqr: Perform inverse quantile regression by solving mixed integer linear programming

#' Perform inverse quantile regression by solving mixed integer linear programming
#' 
#' @param Y data for dependent variable Y
#' @param D data for endogenous variable D
#' @param X data for exogenous variable X
#' @param Z data for instrumental variable Z
#' @param tau tau in quantile regression
#' @param O_pos index for residuals whose sign are fixed as positive
#' @param O_neg index for residuals whose sign are fixed as negative
#' @param eps a small number
#' @param M a large number 
#' @param TimeLimit gurobi parameter
#' @param FeasibilityTol gurobi parameter
#' @param lpsolver "gurobi","cplexapi","lpsolveapi"
#' @return A list containing optimal value, B_D, B_X, B_Z, x and dual variables a 
#' @examples 
#' n=50
#' pD=3
#' sample_data<-chen_lee(n,pD)
#' iqr_milp(Y=sample_data$Y, D=sample_data$D, X = matrix(1, n, 1), Z=sample_data$Z,lpsolver="gurobi") 
#' @description This function is an internal function called inside function "ivqr_fit()"

#####iqr_milp#####
iqr_milp <- function(Y, D, X, Z, tau=0.76, O_pos = NULL, O_neg = NULL, cuts=NULL, eps = 1e-14, M = 10, TimeLimit = 300, FeasibilityTol = 1e-6, lpsolver=NULL){#20191210 cuts is not used?
  
  #### error checking for lpsolver ####
  if (hasArg(lpsolver))
    lpsolver <- tolower(lpsolver)#case-insensitive for user input
  if (is.null(lpsolver)) {
    if (requireNamespace("gurobi", quietly = TRUE)) {
      lpsolver <- "gurobi"
    } else if (requireNamespace("cplexAPI", quietly = TRUE)) {
      lpsolver <- "cplexapi"
    } else if (requireNamespace("lpSolveAPI", quietly = TRUE)) {
      lpsolver <- "lpsolveapi"
    } else {
      stop(
        gsub(
          "\\s+",
          " ",
          "Please install one of the following packages required for
                      estimation:
                      gurobi (version 7.5-1 or later);
                      cplexAPI (version 1.3.3 or later);
                      lpSolveAPI (version 5.5.2.0 or later)."
        ),
        call. = FALSE
      )
    }
  } else {
    if (!lpsolver %in% c("gurobi",
                         "cplexapi",
                         "lpsolveapi")) {
      stop(gsub(
        "\\s+",
        " ",
        paste0(
          "Estimator is incompatible with linear programming
                             package '",
          lpsolver,
          "'. Please install one of the
                             following linear programming packages instead:
                             gurobi (version 7.5-1 or later);
                             cplexAPI (version 1.3.3 or later);
                             lpSolveAPI (version 5.5.2.0 or later)."
        )
      ),
      call. = FALSE)
    } else{
      if (lpsolver == "cplexapi") {
        lpsolver_name = "cplexAPI"
      } else if (lpsolver == "lpsolveapi") {
        lpsolver_name = "lpSolveAPI"
      } else{
        lpsolver_name = lpsolver
      }
      if (!requireNamespace(lpsolver_name, quietly = TRUE)) {
        stop(gsub(
          "\\s+",
          " ",
          paste0(
            "The lpsolver entered ,",
            lpsolver_name,
            " ,is not available. Please
                       enter a lpsolver that has already been installed."
          )
        ),
        call. = FALSE)
      }
    }
  }
  
  p_D = dim(D)[2] ; if(is.null(p_D)){p_D<-0}; if(dim(as.matrix(D))[2]==1){p_D<-1}
  n = dim(X)[1]
  p_X = dim(X)[2]
  p_Z = dim(Z)[2]
  
  ones = rep(1,n)
  
  O = c(O_neg, O_pos)
  #print(paste0("n is ",n))
  #print(paste0("# of fixed residuals is",length(O)))
  if(is.null(O)==F){O <- sort(O)} #put index of pos and neg outlier together
  I = setdiff(1:n, O)
  n_I = length(I)
  
  if(is.null(O)==F){
    a_fix = rep(NA, n); a_fix[O_pos] = 1; a_fix[O_neg] = 0; a_O <- a_fix[O]
    #20191015, add k_fix, l_fix, k_O, l_O
    k_fix = rep(NA, n); k_fix[O_pos] = 1; k_fix[O_neg] = 0; k_O <- k_fix[O]
    l_fix = rep(NA, n); l_fix[O_pos] = 0; l_fix[O_neg] = 1; l_O <- l_fix[O]
  }
  
  
  #B_Z_plus, B_Z_min, B_D, B_X, u, v, a_I, k_I, l_I
  obj_core = c(rep(1, 2*p_Z), rep(0, p_D), rep(0, p_X), rep(0, 2*n), rep(0, n_I), rep(0, n_I+n_I))
  
  #primal feasibility
  #20191015, change to diag(n)
  A_pf = cbind(Z, -Z, D, X, diag(n), -diag(n), matrix(0, n, n_I+n_I+n_I)) 	
  b_pf = Y
  sense_pf = rep("=", n)
  
  #dual feasibility
  A_df_X = cbind(matrix(0, p_X, 2*p_Z+p_D+p_X), matrix(0, p_X, 2*n), t(X[I,,drop=FALSE]), matrix(0, p_X, n_I+n_I)) 
  if(is.null(O)){ b_df_X = (1-tau)*t(X)%*%ones }  
  if(is.null(O)==F){ b_df_X = (1-tau)*t(X)%*%ones - t(X[O,,drop=FALSE])%*%a_O }
  sense_df_X = rep("=", p_X)
  
  A_df_Z = cbind(matrix(0, p_Z, 2*p_Z+p_D+p_X), matrix(0, p_Z, 2*n), t(Z[I,,drop=FALSE]), matrix(0, p_Z, n_I+n_I)) 
  if(is.null(O)){ b_df_Z = (1-tau)*t(Z)%*%ones }
  if(is.null(O)==F){ b_df_Z = (1-tau)*t(Z)%*%ones - t(Z[O,,drop=FALSE])%*%a_O }
  sense_df_Z = rep("=", p_Z)
  
  
  
  #complementary slackness .. do I need the eps? 
  
  #20191015, # of constraints becomes n;  change the mat corresponding to l; change the rhs
  A_vl_I = cbind(matrix(0, n, (2*p_Z + p_D + p_X)), matrix(0,n,n), diag(n), matrix(0, n, n_I), matrix(0, n, n_I), - M*diag(as.numeric(is.element(1:n, I)))[,is.element(1:n, I)]) #v<=lM
  b_vl_I = rep(0 + eps, n)
  if(is.null(O)==F){b_vl_I[O]=M*l_O}
  sense_vl_I = rep("<=", n)   
  
  #20191015, # of constraints becomes n;  change the mat corresponding to k; change the rhs
  A_uk_I = cbind(matrix(0, n, (2*p_Z + p_D + p_X)), diag(n) , matrix(0,n,n), matrix(0, n, n_I), - M*diag(as.numeric(is.element(1:n, I)))[,is.element(1:n, I)], matrix(0, n, n_I)) #u<=kM
  b_uk_I = rep(0 + eps, n)
  if(is.null(O)==F){b_uk_I[O]=M*k_O}
  sense_uk_I = rep("<=", n)  
  
  A_ak_I = cbind(matrix(0, n_I, (2*p_Z + p_D + p_X + 2*n)), diag(n_I), -diag(n_I), matrix(0, n_I, n_I)) #a_i >= k_i 
  b_ak_I = rep(0, n_I) 
  sense_ak_I =  rep(">=", n_I) 
  
  A_al_I = cbind(matrix(0, n_I, (2*p_Z + p_D + p_X +2*n)), diag(n_I), matrix(0, n_I, n_I), diag(n_I)) #a_i <= 1 - l_i 
  b_al_I = rep(1, n_I) 
  sense_al_I =  rep("<=", n_I) 
  
  
  A_core = rbind(A_pf, A_df_X, A_df_Z, A_vl_I, A_uk_I, A_ak_I, A_al_I)
  b_core = c(b_pf, b_df_X, b_df_Z, b_vl_I, b_uk_I, b_ak_I, b_al_I)
  sense_core = c(sense_pf, sense_df_X, sense_df_Z, sense_vl_I, sense_uk_I, sense_ak_I, sense_al_I)
  
  #B_Z_plus, B_Z_min, B_D, B_X, u_Iu, v_Iv, a_I, k_I, l_I
  lb_core = c(rep(0, 2*p_Z), rep(-Inf, p_D+p_X), rep(0, 2*n), rep(0, 3*n_I))
  ub_core = c(rep(Inf, 2*p_Z), rep(Inf, p_D+p_X), rep(Inf, 2*n), rep(1, 3*n_I))
  vtype_core = c(rep("C", 2*p_Z+p_D+p_X+2*n+n_I), rep("B", 2*n_I))
  
  
  model_core = list()
  model_core$obj = obj_core
  model_core$A = A_core
  model_core$rhs = b_core
  model_core$sense = sense_core
  model_core$lb = lb_core
  model_core$ub = ub_core

  params <- list(TimeLimit = TimeLimit, FeasibilityTol = FeasibilityTol, OutputFlag=0)
  
  if (lpsolver == "gurobi") {
    model_core$vtype = vtype_core
    model_core$modelsense = "min"
    out <- gurobi::gurobi(model_core, params) 
  }
  if (lpsolver == "lpsolveapi") {
    out <- runLpSolveAPI(model_core, "min")
  }
  if (lpsolver == "cplexapi") {
    out <- runCplexAPI(model_core, cplexAPI::CPX_MIN)
  }
  
  
  if(out$status == "OPTIMAL"||out$status == "SUBOPTIMAL"){
    objval = out$objval
    x = out$x
    B_Z_plus = x[1:p_Z]
    B_Z_min = x[(p_Z+1):(2*p_Z)]
    if(p_D>0){ B_D = x[(2*p_Z+1):(2*p_Z+p_D)] }; if(p_D==0){ B_D = NULL }
    B_X = x[(2*p_Z+p_D+1):(2*p_Z+p_D+p_X)]
    u_Iu = x[(2*p_Z+p_D+p_X+1):(2*p_Z+p_D+p_X+n)] 
    v_Iv = x[(2*p_Z+p_D+p_X+n+1):(2*p_Z+p_D+p_X+2*n)] 
    a_I = x[(2*p_Z+p_D+p_X+2*n+1):(2*p_Z+p_D+p_X+2*n+n_I)] 
    k_I = x[(2*p_Z+p_D+p_X+2*n+n_I+1):(2*p_Z+p_D+p_X+2*n+2*n_I)] 
    l_I = x[(2*p_Z+p_D+p_X+2*n+2*n_I+1):(2*p_Z+p_D+p_X+2*n+3*n_I)] 
    #put full dual back together
    if(is.null(O)){
      a = a_I
      k = k_I
      l = l_I
    }
    if(is.null(O)==F){ 
      a_fix[I] = a_I; a <- a_fix
      k_fix[I] = k_I; k <- k_fix
      l_fix[I] = l_I; l <- l_fix
    }
    regression_out = list("objval" = objval, "B_Z" =  B_Z_plus-B_Z_min, "B_D" = B_D, "B_X" = B_X, "a" = a, "x" = x, "status" = out$status )
    return(regression_out)
  }else{
    regression_out = list("status" = out$status )
    return(regression_out)
  }
}