R/preprocess_iqr_milp.R

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

#' Preprocess and call iqr_milp
#' 
#' @param Y data for dependent variable Y passed into iqr_milp
#' @param D data for endogenous variable D passed into iqr_milp
#' @param X data for exogenous variable X passed into iqr_milp
#' @param Z data for instrumental variable Z passed into iqr_milp
#' @param tau tau in quantile regression passed into iqr_milp
#' @param eps a small number passed into iqr_milp
#' @param M a large number passed into iqr_milp
#' @param TimeLimit gurobi parameter passed into iqr_milp
#' @param FeasibilityTol gurobi parameter passed into iqr_milp
#' @param r the rate at which we increase bandwidth
#' @param prop_fixed_start the starting paramter for bandwidth
#' @param lpsolver "gurobi","cplexapi","lpsolveapi"
#' @return A list containing the iqr_milp return, the time consumed, the index of residuals whose signs are fixed
#' @examples 
#' n=50
#' pD=3
#' sample_data<-chen_lee(n,pD)
#' preprocess_iqr_milp(Y=sample_data$Y, D=sample_data$D, X = matrix(1, n, 1), Z=sample_data$Z,lpsolver="gurobi") 

preprocess_iqr_milp<- function(Y, D, X, Z, tau=0.5, eps = 1e-14, M = 10, TimeLimit = 300, FeasibilityTol = 1e-6, r=1.1, prop_fixed_start = 0.4,lpsolver=NULL){
  start_clock = Sys.time()      	
  p_D = dim(D)[2] ; if(is.null(p_D)){p_D<-0}
  n = dim(X)[1]
  O_pos = NULL
  O_neg = NULL           
  XX = cbind(D, X)
  qr = quantreg::rq(Y~XX-1)
  fit = qr$fitted.values
  start_width = quantile(abs(qr$res), prop_fixed_start)
  status = "TIME_LIMIT"
  obj     = 0.5
  width = start_width
  while(((status=="TIME_LIMIT")|(obj!=0))){ 
    lb = fit - width
    ub = fit + width
    O_pos = which(Y > ub)
    O_neg = which(Y < lb)
    O = c(O_pos, O_neg)
    if(is.null(O)==F){O = sort(O)}
    I = setdiff(1:n, O)
    n_I = length(I)
    TT  = exp(n_I/200 + p_D/5 + n_I*p_D/1000)*4 #heuristic for time limit
    if(length(O)==0){TT=Inf}
    fit_full = iqr_milp(Y, D, X, Z, tau=0.5, O_pos= O_pos, O_neg= O_neg, eps = eps, M = M, TimeLimit = TT,lpsolver =lpsolver )
    status = fit_full$status
    obj     = fit_full$objval
    if(is.null(obj)){
      obj = 0.5
    }
    width = width*r
  } #end of WHILE loop
  out = list("out" = fit_full, "restart" = log(width/start_width)/log(1/r)) #not gonna work if width is forced to zero             	
  end_clock = Sys.time() 
  return(list("out" = out, "time" = difftime(end_clock,start_clock,units="mins"), "O_pos" = O_pos, "O_neg" = O_neg))                       	
}