# This is the R program for the Monte Carlo of the paper "Subsampling Inference on Quantile Regression Processes" by V. Chernozhukov and I. Fernandez-Val ############################################################################################## # Section I: Preliminary Programs ############################################################################################## ######################## generate data for i-th monte-carlo #################################### gen.data<- function(gamma, nn) { XX<- rnorm(nn); YY<- XX + (1 + gamma*XX)*rnorm(nn); return(cbind(YY,XX)); } ############## computes the inference process for i-th subset of data ############################## estimates<- function(i,ss,X,Y,taus){ x<- X[ss[,i]] y<- Y[ss[,i]] form<-y~x sol<- rq(form,taus); solqr<- sol$coef[2, ]; solols<- lm(y~x)$coef[-1]; inference.process<- solqr-solols return(inference.process) } ############## computes the critical value for i-th monte-carlo ############################## subsample<- function(X, Y, nn=length(Y), b=b, B, taus) { ss<- matrix(sample(nn, b * B, replace = T), b, B); outpt<- lapply (1:B, FUN=estimates, ss, X,Y,taus) return(outpt) } ############################################################################################## # Section II: Monte-Carlo ############################################################################################# library(quantreg); ########### computes an iteration of the Monte Carlo ######################################### montecarlo_iteration <- function(gamma, nn, b, B, taus=taus, filenSM, filenKS, filenSM_U, filenKS_U, ...) { # generate data datt<- gen.data(gamma, nn); Y<- datt[,1]; X<- datt[,2]; form <- Y~X; # compute the inference process for sample size N: sol<- rq(form,taus); solqr<- sol$coef[2, ]; solols<- lm(Y~X)$coef[-1] inference.processN <- solqr-solols # subsampling stage to determine critical value output<- subsample(X, Y, nn, b, B, taus) out<- t(matrix(unlist(output), length(taus), B)); varr<- b*apply(out, 2, var) ; vb<- sqrt( b*(out-matrix(inference.processN, B, length(taus), byrow=T))^2/matrix(varr, B, length(taus), byrow=T)); # CRAMER - VON MISES STATISTIC (CENTERED) KSM<- apply(vb, 1, mean); KS<- sqrt(nn*mean(inference.processN^2/varr) ); # critical<- rq(KSM~1, tau=.95, alpha=.2)$ci[2]; alp <- .05; # Adjusted critical value that accounts for finite sample bias of the subsampling procedure critical<-quantile(KSM,.95) + qnorm(1-alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KSM),quantile(KSM,.95))[1]); # record the KS, critical value, and the decision reject<- ifelse(KS>critical, 1, 0) write.table(t(c(KS, critical, reject)), file=filenSM, sep=" ", append=T, na = NA, row.names = F, col.names = F) # KOLMOGOROV - SMIRNOFF STATISTIC (CENTERED) KSM<- apply(vb, 1, max); KS<- sqrt(nn*max(inference.processN^2/varr) ); # critical<- rq(KSM~1, tau=.95, alpha=.1, method="sparsity", hs=F)$ci[1]; alp <- .05; # Adjusted critical value that accounts for finite sample bias of the subsampling procedure critical<-quantile(KSM,.95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KSM),quantile(KSM,.95))[1]); # record the KS, critical value, and the decision reject<- ifelse(KS>critical, 1, 0) write.table(t(c(KS, critical, reject)), file=filenKS, sep=" ", append=T, na = NA, row.names = F, col.names = F) # CRAMER - VON MISES STATISTIC (UNCENTERED) vb<- sqrt( b*(out)^2/matrix(varr, B, length(taus), byrow=T)); KSM<- apply(vb, 1, mean); KS<- sqrt(nn*mean(inference.processN^2/varr) ); # critical<- rq(KSM~1, tau=.95, alpha=.2)$ci[2]; alp <- .05; # Adjusted critical value that accounts for finite sample bias of the subsampling procedure critical<-quantile(KSM,.95) + qnorm(1-alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KSM),quantile(KSM,.95))[1]); # record the KS, critical value, and the decision reject<- ifelse(KS>critical, 1, 0) write.table(t(c(KS, critical, reject)), file=filenSM_U, sep=" ", append=T, na = NA, row.names = F, col.names = F) # KOLMOGOROV - SMIRNOFF STATISTIC (UNCENTERED) KSM<- apply(vb, 1, max); KS<- sqrt(nn*max(inference.processN^2/varr) ); # critical<- rq(KSM~1, tau=.95, alpha=.1, method="sparsity", hs=F)$ci[1]; alp <- .05; # Adjusted critical value that accounts for finite sample bias of the subsampling procedure critical<-quantile(KSM,.95) + qnorm(alp/2)*sqrt(.95*(1-.95)/B)/unlist(akj(sort(KSM),quantile(KSM,.95))[1]); # record the KS, critical value, and the decision reject<- ifelse(KS>critical, 1, 0) write.table(t(c(KS, critical, reject)), file=filenKS_U, sep=" ", append=T, na = NA, row.names = F, col.names = F) } ################# computes the Monte Carlo calling iteratively the previous function ######################################################## monteSM_KS<- function(R,B, gamma, nn, b, taus=taus) { filenSM<<- paste(c( "H:\\IQR\\Monte_Carlo\\Results\\result.SM.B", B, ".gamma.", gamma, "K", K, "nn", nn, ".R1000.adj5.txt"), collapse=""); filenKS<<- paste(c( "H:\\IQR\\Monte_Carlo\\Results\\result.KS.B", B, ".gamma.", gamma, "K", K, "nn", nn, ".R1000.adj5.txt"), collapse=""); filenSM_U<<- paste(c( "H:\\IQR\\Monte_Carlo\\Results\\result.SM.B", B, ".gamma.", gamma, "K", K, "nn", nn, ".R1000.adj5.unc.txt"), collapse=""); filenKS_U<<- paste(c( "H:\\IQR\\Monte_Carlo\\Results\\result.KS.B", B, ".gamma.", gamma, "K", K, "nn", nn, ".R1000.adj5.unc.txt"), collapse=""); montecarlo<- function(i, R, B, gamma, nn, b, taus=taus) { montecarlo_iteration(gamma, nn, b, B, taus, filenSM=filenSM, filenKS=filenKS, filenSM_U=filenSM_U, filenKS_U=filenKS_U) } ; now <- proc.time()[3] # checkpoint; for(i in 1:R) montecarlo(i, R, B, gamma, nn, b, taus); speed <- (proc.time()[3]-now)/3600; mmmSM<- apply(read.table(filenSM),2, mean)[3] mmmKS<- apply(read.table(filenKS),2, mean)[3] mmmSM_U<- apply(read.table(filenSM_U),2, mean)[3] mmmKS_U<- apply(read.table(filenKS_U),2, mean)[3] return(unlist(speed), mmmSM, mmmKS, mmmSM_U, mmmKS_U); } # Set parameters: grid of quantiles taus <- (1:19)/20; set.seed(8); # Size of the tests: Gamma = 0 K <- 66; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 100, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 400, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 800, b <- ceiling(10+nn^(1/2.01)), taus) # Large sample size properties of the tests K <- 66; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 1000, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 2000, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 5000, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 10000, b <- ceiling(10+nn^(1/2.01)), taus) # Local power of the tests: Gamma = 0.2 K <- 66; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.2, nn<- 100, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.2, nn<- 400, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.2, nn<- 800, b <- ceiling(10+nn^(1/2.01)), taus) # Global power of the tests: Gamma = 0.5 K <- 66; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.5, nn<- 100, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.5, nn<- 400, b <- ceiling(10+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.5, nn<- 800, b <- ceiling(10+nn^(1/2.01)), taus) #################################################################################################### # Size of the tests: Gamma = 0 K <- 77; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 100, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 400, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 800, b <- ceiling(20+nn^(1/2.01)), taus) # Large sample size properties of the tests K <- 77; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 1000, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 2000, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 5000, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.0, nn<- 10000, b <- ceiling(20+nn^(1/2.01)), taus) # Local power of the tests: Gamma = 0.2 K <- 77; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.2, nn<- 100, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.2, nn<- 400, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.2, nn<- 800, b <- ceiling(20+nn^(1/2.01)), taus) # Global power of the tests: Gamma = 0.5 K <- 77; monteSM_KS(R<- 1000, B<- 250, gamma<- 0.5, nn<- 100, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.5, nn<- 400, b <- ceiling(20+nn^(1/2.01)), taus) monteSM_KS(R<- 1000, B<- 250, gamma<- 0.5, nn<- 800, b <- ceiling(20+nn^(1/2.01)), taus) R<- 250; B<- 100; gamma<- 0.2; K<-5; nn<- 100; Stat<-"KS" R<- 1000; B<- 250; gamma<- 0; K<- 5; nn<- 400; Stat<-"KS" mmm<-read.table(paste(c( "c:\\aaa\\comput\\BQR\\result", Stat, ".gamma.", gamma, "K", K, "nn", nn, ".txt"), collapse="")); apply(mmm,2, mean) mmm<-read.table(paste(c( "c:\\aaa\\comput\\BQR\\result", Stat, "B", B, ".gamma.", gamma, "K", K, "nn", nn, ".txt"), collapse="")); apply(mmm,2, mean)