# this section implments the estimation and inference results 

source("c:\\aaa\\comput\\rq.s"); dyn.load("c:\\aaa\\comput\\rq.obj");

options(object.size=250000000)


#




# section 1:  General Results


Y<- data[1527:2527,5]
X<- data[1527:2527,2:4]

y0<- Y; x0<- X;

formula <- y0 ~ x0;  p <- dim(X)[2]+1



bandwidth.rq<-  function (p, n, hs = TRUE, alpha = 0.05) 
{
    PI <- pi
    x0 <- qnorm(p)
    f0 <- (1/sqrt(2 * PI)) * exp(-((x0^2/2)))
    if (hs == TRUE) 
        n^(-1/3) * qnorm(1 - alpha/2)^(2/3) * ((1.5 * f0^2)/(2 * 
            x0^2 + 1))^(1/3)
    else n^-0.2 * ((4.5 * f0^4)/(2 * x0^2 + 1)^2)^0.2
}

tau<- .01

   h <- bandwidth.rq(p=tau, n=1000, hs = T)
        if (tau + h > 1) 
            stop("tau + h > 1:  error in summary.rq")
        if (tau - h < 0) 
            stop("tau - h < 0:  error in summary.rq")
        uhat <- Y
        h <- (qnorm(tau + h) - qnorm(tau - h))*min(sqrt(var(uhat)), ( quantile(uhat,.75)- quantile(uhat, .25))/1.34 )





taus <- (1:199)/200; fit <- array(0,c(p,4,length(taus)))

for(i in 31:length(taus)){
       print(taus[i])
	robj<- rq(x0, y0,tau=taus[i], method="score", alpha=.1, hs=F);
       fit[,1,i] <- robj$coef ;
      	fit[,2,i] <- robj$ci[1,];
	fit[,3,i] <- robj$ci[2,];
        }


tausl <- (1:30)/200; 

for(i in 1:length(tausl)){
       print(tausl[i])
	robj<- rq(x0, y0,tau=tausl[i], method="score", alpha=.1, hs=F);
       fit[,1,i] <- robj$coef ;
       ci<-  rq.infer.extreme(X,Y, tau=tausl[i], psub=.25, space=20, B=100);
       fit[,2,i] <- ci$ciL.e;
	fit[,3,i] <- ci$ciU.e;
        }



tausu <- (170:199)/200;
for(i in 1:30){
       print(tausu[i]); 
      	robj<- rq(x0, -y0, tau=(1- tausu[i]), method="score", alpha=.1, hs=F);
       fit[,1,169+i] <- -robj$coef ;
       ci<-  rq.infer.extreme(x0,-y0, tau=(1- tausu[i]), psub=.25, space=20, B=100);
       fit[,2,169+i] <- -ci$ciL.e;
	fit[,3,169+i] <- -ci$ciU.e;
        }





postscript("c:\\aaa\\\comput\\cec\\result1.ps",horizontal=F,font=5,pointsize=11,width=4.0,height=7.0)


blab_c("Intercept", "Spot Oil Return", "Market Return", "Own Return")



par(mfrow=c(2,2))

for(i in 1:p){
        b<-   fit[i,1,];
        b.p<- fit[i,3,]; 
        b.m<- fit[i,2,]; 
        plot( c(0,1 ), range(c(b, b.m,b.p)), xlim =c(0, 1), type="n",xlab="Quantile Index", ylab="Coefficient", main=blab[i]);
        polygon(c(taus,rev(taus)),c(b.m,rev(b.p)),density=-5, border=T, col=2);
        lines(taus,b);
	 
        }
dev.off();



 
tausl <- (1:30)/200;  fitLK <- array(0,c(p,4,length(tausl))) ; fitLB <- array(0,c(p,4,length(tausl)))
tausl <- (1:30)/200;  fitUK <- array(0,c(p,4,length(tausl))) ; fitUB <- array(0,c(p,4,length(tausl)))


for(i in 1:length(tausl)){
       print(tausl[i])
	robj<- rq(x0, y0,tau=tausl[i], method="score", alpha=.1, hs=F);
       fitLK[,1,i] <- robj$coef ;   	fitLB[,1,i] <- robj$coef ;
       ci<-  rq.infer.central(X,Y, tau=tausl[i], ban=.015, B=250) ;
       fitLK[,2,i] <- ci$ciL.ck;    fitLB[,2,i] <- ci$ciL.cb;
 	fitLK[,3,i] <- ci$ciU.ck;    fitLB[,3,i] <- ci$ciU.cb;  }


for(i in 1:length(tausl)){
       print(tausl[i])
	robj<- rq(x0, -y0, tau=tausl[i], method="score", alpha=.1, hs=F);
       fitUK[,1,i] <- -robj$coef ;   fitUB[,1,i] <- robj$coef ;
       ci<-  rq.infer.central(x0,-y0, tau=tausl[i], ban=.015, B=250) ;
       fitUK[,2,i] <- -ci$ciL.ck;    fitUB[,2,i] <- -ci$ciL.cb;
	fitUK[,3,i] <- -ci$ciU.ck;    fitUB[,3,i] <- -ci$ciU.cb;
        }





postscript("c:\\aaa\\\comput\\cec\\result2.ps",horizontal=F,font=5,pointsize=11,width=4.0,height=7.0)


par(mfrow=c(2,2))

for(i in 1:p){
        b<-   fitLK[i,1,];
        bK.p<- fitLK[i,3,]; 
        bK.m<- fitLK[i,2,]; 
	 be.p<- fit[i,3,1:30]; 
        be.m<- fit[i,2,1:30]; 
	 bB.p<- fitLB[i,3,]; 
        bB.m<- fitLB[i,2,]; 
        plot( c(0,.16 ), range(c(b, be.m,be.p)), xlim =c(0, .16), type="n",xlab="Quantile Index", ylab="Coefficient", main=blab[i]);
        lines(tausl,b);
        points(tausl,b);
	 lines(tausl, bK.p, lty=2);
	 lines(tausl, bK.m, lty=2);
	 lines(tausl,smooth(be.p), lty=1);
	 lines(tausl,smooth(be.m), lty=1);

        }
dev.off();






postscript("c:\\aaa\\\comput\\cec\\result3.ps",horizontal=F,font=5,pointsize=11,width=4.0,height=7.0)




par(mfrow=c(2,2))

for(i in 1:p){
        b<-   fitLK[i,1,];
        bK.p<- fitLK[i,3,]; 
        bK.m<- fitLK[i,2,]; 
	 be.p<- fit[i,3,1:30]; 
        be.m<- fit[i,2,1:30]; 
	 bB.p<- fitLB[i,3,]; 
        bB.m<- fitLB[i,2,]; 
        plot( c(0,.16 ), range(c(b, be.m,be.p)), xlim =c(0, .16), type="n",xlab="Quantile Index", ylab="Coefficient", main=blab[i]);
        lines(tausl,b);
        points(tausl,b);
	 lines(tausl,bB.p, lty=2);
	 lines(tausl,bB.m, lty=2);
	 lines(tausl,smooth(be.p), lty=1);
	 lines(tausl,smooth(be.m), lty=1);

        }
dev.off();






postscript("c:\\aaa\\\comput\\cec\\result4.ps",horizontal=F,font=5,pointsize=11,width=4.0,height=7.0)


par(mfrow=c(2,2))

for(i in 1:p){
        b<-   fitUK[i,1,];
        bK.p<- fitUK[i,3,]; 
        bK.m<- fitUK[i,2,]; 
	 be.p<- fit[i,3,170:199]; 
        be.m<- fit[i,2,170:199]; 
        plot( c(.84, 1 ), range(c(b,bK.p, bK.m, be.m,be.p)), xlim =c(.84, 1), type="n",xlab="Quantile Index", ylab="Coefficient", main=blab[i]);
        lines(1- tausl,b);
        points(1-tausl,b);
	 lines(1-tausl,bK.p, lty=2);
	 lines(1-tausl,bK.m, lty=2);
	 lines(.845+ tausl,smooth(be.p), lty=1);
	 lines(.845+ tausl,smooth(be.m), lty=1);

        }
dev.off();






postscript("c:\\aaa\\\comput\\cec\\result5.ps",horizontal=F,font=5,pointsize=11,width=4.0,height=7.0)


par(mfrow=c(2,2))

for(i in 1:p){
        b<-   fitUK[i,1,];
        bK.p<- fitUB[i,3,]; 
        bK.m<- fitUB[i,2,]; 
	 be.p<- fit[i,3,170:199]; 
        be.m<- fit[i,2,170:199]; 
        plot( c(.84, 1 ), range(c(b,bK.p, bK.m, be.m,be.p)), xlim =c(.84, 1), type="n",xlab="Quantile Index", ylab="Coefficient", main=blab[i]);
        lines(1- tausl,b);
        points(1-tausl,b);
	 lines(1-tausl,bK.p, lty=2);
	 lines(1-tausl,bK.m, lty=2);
	 lines(.845+ tausl,smooth(be.p), lty=1);
	 lines(.845+ tausl,smooth(be.m), lty=1);

        }
dev.off();