R/ivqr_fit.R

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

#' IVQR Inference
#'
#' @param 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 M a large number
#' @param homoskedastic
#' @param kernel
#' @param lpsolver "gurobi","cplexapi","lpsolveapi"
#' @return The lower bound and upper bound for B_D and B_X
#' @examples 
#' n=50
#' pD=3
#' sample_data<-chen_lee(n,pD)
#' Y=sample_data$Y
#' D=sample_data$D
#' Z=sample_data$Z
#' X = matrix(1, n, 1)
#' ivqr_fit (Y=Y,D=D,X=X,Z=Z,tau = 0.5,M = 10,homoskedastic = FALSE,kernel = "Powell", lpsolver="gurobi")

ivqr_fit <- function(Y,
                     D,
                     X,
                     Z,
                     tau = 0.5,
                     M = 10,
                     homoskedastic = FALSE,
                     kernel = "Powell",
                     lpsolver = NULL) {

  #### TBA: error checking for all the input parameters ####
  
  iqr_milp_fit = iqr_milp(Y, D, X, Z, tau = tau, lpsolver = lpsolver)
  
  result_list <- list()
  XZ = cbind(X, Z)
  PI_XZ = XZ %*% solve(t(XZ) %*% XZ) %*% t(XZ)  #PI = X%*%solve(t(X)%*%X)%*%t(X)
  
  Phi_CH = PI_XZ %*% D
  
  B_D_hat = iqr_milp_fit$B_D
  for (J in 1:dim(D)[2]) {
    result_list[[J]] <-
      ivqr_confint(
        test_stat = test_stat_freq_D,
        j = J,
        B = Phi_CH[, J, drop = FALSE],
        Phi = Phi_CH[, -J, drop = FALSE],
        Beta_j_hat = B_D_hat[J],
        Srate = 1 / 2,
        alpha = 0.1,
        Y,
        D,
        X,
        Z,
        tau = 0.5,
        M = 10,
        homoskedastic = homoskedastic,
        kernel = kernel,
        lpsolver = lpsolver
      )
    print(paste("J", J))
  }
  
  B_X_hat = iqr_milp_fit$B_X
  for (K in 1:dim(X)[2]) {
    result_list[[K + dim(D)[2]]] <-
      ivqr_confint(
        test_stat = test_stat_freq_X,
        j = K,
        B = NULL,
        Phi = Phi_CH,
        Beta_j_hat = B_X_hat[K],
        Srate = 1 / 2,
        alpha = 0.1,
        Y,
        D,
        X,
        Z,
        tau = 0.5,
        M = 10,
        homoskedastic = homoskedastic,
        kernel = kernel,
        lpsolver = lpsolver
      )
    print(paste("K", K))
  }
  
  z = list()
  z$coefficients <- do.call(rbind, result_list)
  #dimnames(z$coefficients)[[1]]=paste0("B_D_",1:dim(D)[2])
  dimnames(z$coefficients)[[1]] = c(paste0("B_D_", 1:dim(D)[2]), paste0("B_X_", 1:dim(X)[2]))
  dimnames(z$coefficients)[[2]] = c("coefficients", "lower bd", "upper bd")
  z$tau = tau
  z$kernel = kernel
  z$homoskedastic = homoskedastic
  z$call = match.call()
  class(z) <- "ivqr_fit"
  z
}
#' Print the IVQR inference results
#'
#' @param x what ivqr_fit() returns
#' @return the point estimation for the exogenous varialbes and endogenous variables
#' @examples 
#' n=50
#' pD=3
#' sample_data<-chen_lee(n,pD)
#' Y=sample_data$Y
#' D=sample_data$D
#' Z=sample_data$Z
#' X = matrix(1, n, 1)
#' print(ivqr_fit(Y=Y,D=D,X=X,Z=Z,tau = 0.5,M = 10,homoskedastic = FALSE,kernel = "Powell", lpsolver="gurobi"))

print.ivqr_fit <- function (x, ...)
{
  if (!is.null(cl <- x$call)) {
    cat("Call:\n")
    dput(cl)
  }
  
  coef <- coef(x)
  cat("\nCoefficients:\n")
  print(format(round(coef[, 1], digits = 5)), quote = FALSE, ...)
  if (!is.null(attr(x, "na.message")))
    cat(attr(x, "na.message"), "\n")
  invisible(x)
}
#' Print a summary of the IVQR inference results
#'
#' @param x what ivqr_fit() returns
#' @return the point estimation and its confidence interval for the exogenous varialbes and endogenous variables
#' @examples 
#' n=50
#' pD=3
#' sample_data<-chen_lee(n,pD)
#' Y=sample_data$Y
#' D=sample_data$D
#' Z=sample_data$Z
#' X = matrix(1, n, 1)
#' summary(ivqr_fit(Y=Y,D=D,X=X,Z=Z,tau = 0.5,M = 10,homoskedastic = FALSE,kernel = "Powell", lpsolver="gurobi"))

summary.ivqr_fit <- function (x, ...)
{
  if (!is.null(cl <- x$call)) {
    cat("Call:\n")
    dput(cl)
  }
  tau <- x$tau
  cat("\ntau: ")
  print(format(round(tau, digits = 5)), quote = FALSE, ...)
  
  kernel <- x$kernel
  cat("\nkernel: ")
  print(kernel, quote = FALSE, ...)
  
  homoskedastic <- x$homoskedastic
  cat("\nhomoskedastic: ")
  print(homoskedastic, quote = FALSE, ...)
  
  coef <- coef(x)
  #names(coef) <- c("coefficients", "lower bd", "upper bd")
  cat("\nCoefficients:\n")
  print(format(round(coef, digits = 5)), quote = FALSE, ...)
  if (!is.null(attr(x, "na.message")))
    cat(attr(x, "na.message"), "\n")
  invisible(x)
}