#	This is is a collection of very simple functions to perform extremal inference
#	and construct bias-corrected estimates 	by subsampling self-normalized regression quantile statistics.
#	Theory of the method is given in the papers by V.Chernozhukov
#	Project: "Conditional Extremes and Near-Extremes" e-mail: vchern@mit.edu
#	In order to use this one will need a basic quanile regression routine e.g one for windows
#	source("c:\\aaa\\comput\\rq.s"); dyn.load("c:\\aaa\\comput\\rq.obj");
# 	This is the function that does the inference and constructs bias-corrected estimates. Usage:

ci<- rq.infer.extreme(X,Y, tau=.005, psub=.3, space=20, B=100);
ci.c<- rq.infer.central(X,Y, tau=.2, ban=.0003)
#tsplot( ci$ciU.e, ci$ciL.e, ci$BC.e);


rq.infer.extreme<- function(X,Y, tau=.1, psub=.2, space=20, B=200) {
	# X is  n x k regressor matrix
	# Y is  n x 1 vector 
	# Space equals to "tau * (m-1) * n" that should greater than dim(X) ; Inference results should not depend much on this parameter
	# psub determines the subsample size using formula  b = psub * n
	# B is the number of bootstrap draws
	n<- length(Y); 	b<- floor(n/4); 	
	tau.e<- tau;  	
	m<- (tau.e*n+space)/(tau.e*n);   	p <- dim(X)[2] +1
	Res<- null();
	betat<- rq(X,Y, tau=tau.e)$coef ;   
	betatm<- rq(X,Y, tau=m*tau.e)$coef ; 
	muX<- c(1,apply(X,2, mean)) ;  
	An<-  (m-1)*tau*sqrt(n/(tau*(1-tau)) )/(muX%*%(  betatm - betat   ) )
# choice of self-normaziation An is proprtional to that in my papers; done in this way to produce stable results
	tau.b.e<- min(tau*n/b,.3);
# this is the quantile index for computing the order statsitics in the subsamples. It is truncated not to exceed .3
	betatb.e<- rq(X,Y, tau=tau.b.e)$coef ;   
	B<- B;
# set.seed(1);
	 ss<- matrix(sample(1:n, b * B, replace = T), b, B);
	 for(i in 1:B) { 		
		Xs<- X[ss[,i],]; Ys<- Y[ss[,i]];
 		m.b.e<- (tau.b.e*b+space)/(tau.b.e*b);  
		sub.betatb.e<- rq(Xs,Ys, tau=tau.b.e)$coef ;  sub.betatbm.e<- rq(Xs,Ys, tau=m.b.e*tau.b.e)$coef ;
		sub.betatb.c<- rq(Xs,Ys, tau)$coef 
 		Ab.e<- (m.b.e-1)*tau.b.e*sqrt(b/(tau.b.e* (1-tau.b.e)) )/ ( -muX%*%(   sub.betatb.e -   sub.betatbm.e   ));
 		Sb.e<-  Ab.e * (sub.betatb.e -  betatb.e);    
		Sb<- c(Sb.e)
		Res<- rbind(Res, Sb)
	 }
 	ciL.e<- betat - apply(Res[,1:p], 2, quantile, .94)/An ; 
	ciU.e<- betat - apply(Res[,1:p], 2, quantile, .06)/An ;
# Median-Bias Correction:
	BC.e<-  betat - apply(Res[,1:p], 2, quantile, .5)/An ;
	if (tau>.5) {return(-1*c(ciL.e, ciU.e, BC.e, betat)); } else { return(ciL.e, ciU.e, BC.e, betat);}

 }




rq.infer.central<- function(X,Y, tau=.1, ban=.005, B=10) {
	# X is  n x k regressor matrix
	# Y is  n x 1 vector 
	# B is the number of bootstrap draws
	n<- length(Y);   p <- dim(X)[2] +1
	Coef<- null();
	betat<- rq(X,Y, tau=tau)$coef ;  
	Xt<- cbind(rep(1, n), X);  Xt<- as.matrix(Xt)
	J<-  (1/n*1/ban )*t(Xt)%*%diag(dnorm( (Y - t(betat)%*%t(Xt))/ban , mean=0, sd=1) )%*%Xt ;
	Q<- (tau)*(1-tau)*t(Xt)%*%Xt/n ;
	serror<- sqrt(diag(solve(J)%*%Q%*%solve(J))/n);
	ciU.ck<- betat + serror*1.69 ; 
	ciL.ck<- betat - serror*1.69 ;
 	return(ciL.ck, ciU.ck, betat);
 	}