tests/testthat/test20190828.R

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

library(ivqr)
library(quantreg); 
library(MASS)

#########Simulation
n=150 #need to extend time limit argument in iqr_milp function in  for larger n's

Z1 = rnorm(n); Z2 = rnorm(n); Z3 = rnorm(n)

Sigma = matrix(c(1, 0.4, 0.6, -0.2, 0.4, 1, 0, 0, 0.6, 0, 1, 0, -0.2, 0, 0, 1), 4, 4, byrow=TRUE)
EpsV1V2V3 = mvrnorm(n, mu = c(0,0,0,0), Sigma = 0.25*Sigma)

D1D2D3 = matrix(NA, n, 3)
for(i in 1:n){
  D1D2D3[i,1] = pnorm(Z1[i] + EpsV1V2V3[i,2]) #V2 correlate with V1, indep of the others
  D1D2D3[i,2] = 2*pnorm(Z2[i] + EpsV1V2V3[i,3]) #V3 correlate with V1, indep of the others
  D1D2D3[i,3] = 1.5*pnorm(Z3[i] + EpsV1V2V3[i,4]) #V4 correlate with V1, indep of the others
}

Y = 1 + D1D2D3[,1] + D1D2D3[,2] +  D1D2D3[,3] + (0.5 + D1D2D3[,1] + 0.25*D1D2D3[,2] +  0.15*D1D2D3[,3])*EpsV1V2V3[,1]
Y = matrix(Y, n, 1)

X = matrix(1, n, 1) #cbind(rep(1,n), matrix(rnorm(3*n), n, 3)) 
Z = cbind(Z1, Z2, Z3)

D = D1D2D3

rq(Y ~ D + X -1)
result_before<-rq(Y-D%*%rep(1,3) ~ X + Z-1) #for true values of beta_D the instruments have quite small estimated regression coefficients
test_that("not equal to 0 before we find the right D",expect_false(isTRUE(all.equal(unname(result_before$coefficients[2:4]),rep(0,3)))))


#for tau=0.5, true value of beta_D is 1 for all entries 

iqr_milp_fit = iqr_milp(Y, D, X, Z,lpsolver="gurobi") 

iqr_milp_fit$B_D #exact IQR coefficients
iqr_milp_fit$B_Z #Instrumental variables coefficients, zero when just identified
iqr_milp_fit$a #access dual variables directly, they are used in the paper for regression rankscore inference 

result_after<-rq(Y-D%*%iqr_milp_fit$B_D ~ X + Z-1) #compare with the case where B_D=rep(1,3), the coef for Z are close to 0 now
#all.equal(result_after$coefficients[2:4],rep(0,3),check.names=FALSE)
test_that("equal to 0 after we find the right D",expect_equal(unname(result_after$coefficients[2:4]),rep(0,3)))