#The MC is done in S-Plus. Graphics is done in R.
#Y is generated according to a Weibull-error model.
#X is generated by sampling with replacement from var_data.txt 




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

################### Section I: Preliminary Programs ###############################

######################## generate data for i-th monte-carlo ####################################



gen.data.W<- function(nn, DF) {
	ind<- sample(1:1001, nn, replace=T);
    	XX<- as.matrix(Xt[ind,]);
	YY<-  XX%*%rep(1, 4) + .07*rweibull(nn, DF);
	return(cbind(YY,XX) ); }

########################### Monte-Crlo Function #####



filenD<-  "c:\\aaa\\comput\\CEC\\nn500.CASED.HS.run1.txt"
filen.biasD<-  "c:\\aaa\\comput\\CEC\\nn500.CASED.HS.BiasAnalysis.run1.txt"
filenE<-  "c:\\aaa\\comput\\CEC\\nn500.CASEE.HS.run1.txt"
filen.biasE<-  "c:\\aaa\\comput\\CEC\\nn500.CASEE.HS.BiasAnalysis.run1.txt"
filenF<-  "c:\\aaa\\comput\\CEC\\nn500.CASEF.HS.run1.txt"
filen.biasF<-  "c:\\aaa\\comput\\CEC\\nn500.CASEF.HS.BiasAnalysis.run1.txt"


nn<- 500;
Xt<- cbind(rep(1, 1001),X);


monte.cov.W<- function(filen, filen.bias, DF, ...) {

# generate data ; select badnwidth

datt<- gen.data.W(nn, DF=DF);  YY<- datt[,1]; XX<- datt[,-c(1:2)];
dd<- rq(XX,YY, tau=.5)$res
ban<-  ((quantile(dd,.75) - quantile(dd, .25))/1.34)*(500)^{-1/3}  # A Hall-Sheather Rule

# extreme inference

extr.005<- rq.infer.extreme(XX,YY, tau=.005, psub=.2, space=10, B=200);
extr.020<- rq.infer.extreme(XX,YY, tau=.020, psub=.2, space=10, B=200);
extr.050<- rq.infer.extreme(XX,YY, tau=.050, psub=.2, space=10, B=200);
extr.100<- rq.infer.extreme(XX,YY, tau=.100, psub=.2, space=10, B=200);

mise.005<- c(  extr.005$BC.e[1] - 1- .07*qweibull(.005,DF),  extr.005$BC.e[3] - 1);
mise.020<- c(  extr.020$BC.e[1] - 1- .07*qweibull(.020,DF),  extr.020$BC.e[3] - 1);
mise.050<- c(  extr.050$BC.e[1] - 1- .07*qweibull(.010,DF),  extr.050$BC.e[3] - 1);
mise.100<- c(  extr.100$BC.e[1] - 1- .07*qweibull(.100,DF),  extr.100$BC.e[3] - 1);
 

cove.005<-  c( I( extr.005$ciL.e[1] < 1+ .07*qweibull(.005,DF) )* I ( 1+ .07*qweibull(.005, DF)< extr.005$ciU.e[1]), I( extr.005$ciL.e[3] < 1)* I (1< extr.005$ciU.e[3]))
cove.020<-  c( I( extr.020$ciL.e[1] < 1+ .07*qweibull(.020,DF) )* I ( 1+ .07*qweibull(.020, DF)< extr.020$ciU.e[1]), I( extr.020$ciL.e[3] < 1)* I (1< extr.020$ciU.e[3]))
cove.050<-  c( I( extr.050$ciL.e[1] < 1+ .07*qweibull(.050,DF) )* I ( 1+ .07*qweibull(.050, DF)< extr.050$ciU.e[1]), I( extr.050$ciL.e[3] < 1)* I ( 1< extr.050$ciU.e[3]))
cove.100<-  c( I( extr.100$ciL.e[1] < 1+ .07*qweibull(.100,DF) )* I ( 1+ .07*qweibull(.100, DF)< extr.100$ciU.e[1]), I( extr.100$ciL.e[3] < 1)* I ( 1< extr.100$ciU.e[3]))

centr.005<- rq.infer.central(XX,YY, tau=.005, ban=ban)
centr.020<- rq.infer.central(XX,YY, tau=.020, ban=ban)
centr.050<- rq.infer.central(XX,YY, tau=.050, ban=ban)
centr.100<- rq.infer.central(XX,YY, tau=.100, ban=ban)

covc.005<-  c( I( centr.005$ciL.ck[1] < 1+ .07*qweibull(.005,DF) )* I ( 1+ .07*qweibull(.005,DF)< centr.005$ciU.ck[1]), I( centr.005$ciL.ck[3] < 1)* I ( 1< centr.005$ciU.ck[3]))
covc.020<-  c( I( centr.020$ciL.ck[1] < 1+ .07*qweibull(.020,DF) )* I ( 1+ .07*qweibull(.020,DF)< centr.020$ciU.ck[1]), I( centr.020$ciL.ck[3] < 1)* I ( 1< centr.020$ciU.ck[3]))
covc.050<-  c( I( centr.050$ciL.ck[1] < 1+ .07*qweibull(.050,DF) )* I ( 1+ .07*qweibull(.050,DF)< centr.050$ciU.ck[1]), I( centr.050$ciL.ck[3] < 1)* I ( 1< centr.050$ciU.ck[3]))
covc.100<-  c( I( centr.100$ciL.ck[1] < 1+ .07*qweibull(.100,DF) )* I ( 1+ .07*qweibull(.100,DF)< centr.100$ciU.ck[1] ), I( centr.100$ciL.ck[3] <1)* I ( 1< centr.100$ciU.ck[3]))

misc.005<- c(  centr.005$betat[1] - 1- .07*qweibull(.005,DF) ,  centr.005$betat[3] - 1 );
misc.020<- c(  centr.020$betat[1] - 1- .07*qweibull(.020,DF) ,  centr.020$betat[3] - 1 );
misc.050<- c(  centr.050$betat[1] - 1- .07*qweibull(.010,DF) ,  centr.050$betat[3] - 1 );
misc.100<- c(  centr.100$betat[1] - 1- .07*qweibull(.100,DF) ,  centr.100$betat[3] - 1 );
 

cov.extr<-  c( cove.005, cove.020, cove.050, cove.100) ;
miss.extr<- c( mise.005, mise.020, mise.050, mise.100) ;
cov.centr<- c( covc.005, covc.020, covc.050, covc.100) ;
miss.centr<- c( misc.005, misc.020, misc.050, misc.100) ;


write.table(t(c(cov.extr, cov.centr)), file=filen, sep="  ", append=T,  dimnames.write = F, na = NA)
write.table(t(c(miss.extr, miss.centr)), file=filen.bias, sep="  ", append=T,  dimnames.write = F, na = NA)

}


now<- proc.time(); for(i in 1:100) {monte.cov.W(filenD, filen.biasD, DF=.5 )} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:100) {monte.cov.W(filenE, filen.biasE, DF=1  )} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:100) {monte.cov.W(filenF, filen.biasF, DF=3)} ;   speed <- (proc.time()-now)/3600;


now<- proc.time(); for(i in 1:100) {monte.cov.W(filenD, filen.biasD, DF=.5 )} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:100) {monte.cov.W(filenE, filen.biasE, DF=1  )} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:100) {monte.cov.W(filenF, filen.biasF, DF=3)} ;   speed <- (proc.time()-now)/3600;


now<- proc.time(); for(i in 1:100) {monte.cov.W(filenD, filen.biasD, DF=.5 )} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:100) {monte.cov.W(filenE, filen.biasE, DF=1  )} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:100) {monte.cov.W(filenF, filen.biasF, DF=3)} ;   speed <- (proc.time()-now)/3600;




# This part is DONE in R:

# NOTE YET !!!!!!!!!!!!!!!!


volD<- apply(read.table(filenD),2, mean)
volE<- apply(read.table(filenE),2, mean)
volF<- apply(read.table(filenF),2, mean)

taus<- c(.005, .02, .05, .1)


postscript("c:\\aaa\\\comput\\cec\\MCresultDEF.ps",horizontal=F,pointsize=9,width=7.0,height=7.0)

par(mfrow=c(3,2), lab=c(10,10,10))

plot( range(taus), range(0,1), type="n", ylim=c(0,1), xlim=c(0, .11), ylab="Coverage", xlab="Tau ", main="D. Weibull(.5)-disturbance: Coverage for Intercept", lwd=2)
lines( taus, volD[seq(1,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volD[seq(9,16, by=2)], type="l", lty=4, lwd=2, col=2)
abline(h=.9, col=1, lwd=2, lty=2)
typ.names <- c("Extremal ", "Central", "Nominal 90%")
legend(.05, .4, legend = typ.names, lty = c(1,4,2), col=c(4,2,1), bty="no", lwd=2)

plot( range(taus), range(0,1), type="n", ylim=c(0,1), xlim=c(0, .11), ylab="Coverage", xlab="Tau ", main="D.  Weibull(.5)-disturbance: Coverage for Slopes", lwd=2)
lines( taus, volD[seq(2,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volD[seq(10,16, by=2)], type="l", lty=4, lwd=2, col=2)
abline(h=.9, col=1, lwd=2, lty=2)
typ.names <- c("Extremal ", "Central", "Nominal 90%")
legend(.05, .4, legend = typ.names, lty = c(1,4,2), col=c(4,2,1), bty="no", lwd=2)


plot( range(taus), range(0,1), type="n", ylim=c(0,1), xlim=c(0, .11), ylab="Coverage", xlab="Tau ", main="E.  Weibull(1)-disturbance: Coverage for Intercept", lwd=2)
lines( taus, volE[seq(1,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volE[seq(9,16, by=2)], type="l", lty=4, lwd=2, col=2)
abline(h=.9, col=1, lwd=2, lty=2)
typ.names <- c("Extremal ", "Central", "Nominal 90%")
legend(.05, .4, legend = typ.names, lty = c(1,4,2), col=c(4,2,1), bty="no", lwd=2)



plot( range(taus), range(0,1), type="n", ylim=c(0,1), xlim=c(0, .11), ylab="Coverage", xlab="Tau ", main="E.  Weibull(1)-disturbance: Coverage for Slopes", lwd=2)
lines( taus, volE[seq(2,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volE[seq(10,16, by=2)], type="l", lty=4, lwd=2, col=2)
abline(h=.9, col=1, lwd=2, lty=2)
typ.names <- c("Extremal ", "Central", "Nominal 90%")
legend(.05, .4, legend = typ.names, lty = c(1,4,2), col=c(4,2,1), bty="no", lwd=2)



plot( range(taus), range(0,1), type="n", ylim=c(0,1), xlim=c(0, .11), ylab="Coverage", xlab="Tau ", main="F.  Weibull(4)-disturbance: Coverage for Intercept", lwd=2)
lines( taus, volF[seq(1,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volF[seq(9,16, by=2)], type="l", lty=4, lwd=2, col=2)
abline(h=.9, col=1, lwd=2, lty=2)
typ.names <- c("Extremal ", "Central", "Nominal 90%")
legend(.05, .4, legend = typ.names, lty = c(1,4,2), col=c(4,2,1), bty="no", lwd=2)


plot( range(taus), range(0,1), type="n", ylim=c(0,1), xlim=c(0, .11), ylab="Coverage", xlab="Tau ", main="F.  Weibull(4)-disturbance: Coverage for Slopes", lwd=2)
lines( taus, volF[seq(2,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volF[seq(10,16, by=2)], type="l", lty=4, lwd=2, col=2)
abline(h=.9, col=1, lwd=2, lty=2)
typ.names <- c("Extremal ", "Central", "Nominal 90%")
legend(.05, .4, legend = typ.names, lty = c(1,4,2), col=c(4,2,1), bty="no", lwd=2)



dev.off();








medad<- function(x){ median(abs(x))} 
absmed<- function(x){ abs(median(x))} 

volD<- apply(read.table(filen.biasD),2, medad)
volE<- apply(read.table(filen.biasE),2, medad)
volF<- apply(read.table(filen.biasF),2, medad)


#plot(abs(volD[1:8])/ abs(volD[9:16]))



postscript("c:\\aaa\\\comput\\cec\\MC.Forecast.resultDEF.ps",horizontal=F,pointsize=9,width=7.0,height=7.0)

par(mfrow=c(3,2), lab=c(10,10,10))

xlab<- "Tau"
ylab<- "Median Abs Deviation"
xlim<- c(0.005, .1)
ylim<- c(.5,1.5)

plot( range(.005,.1), range(.5,1.5), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="A. Weibull(.5)-disturbance: Estimating Intercept", lwd=2)
lines( taus, c(volD[seq(1,8, by=2)]/volD[seq(9,16, by=2)]), type="l", lty=2, lwd=2, col=3)
typ.names <- c("Ratio of MAD of Bias-Corrected QR \n to MAD of QR")
abline(h=1, col=1, lwd=2, lty=1)
legend(.025, .7, legend = typ.names, lty = c(2), col=c(3), bty="no", lwd=2)


plot( range(.005,.1), range(0,60), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="A. Weibull(.5)-disturbance: Estimating Slopes", lwd=2)
lines( taus, volD[seq(2,8, by=2)]/volD[seq(10,16, by=2)], type="l", lty=2, lwd=2, col=3)
typ.names <- c("Ratio of MAD of Bias-Corrected QR \n to MAD of QR")
abline(h=1, col=1, lwd=2, lty=1)
legend(.025, .7, legend = typ.names, lty = c(2), col=c(3), bty="no", lwd=2)


plot( range(.005,.1), range(0,60), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="B. Weibull(1)-disturbance: Estimating Intercept", lwd=2)
lines( taus, volE[seq(1,8, by=2)]/volE[seq(9,16, by=2)], type="l", lty=2, lwd=2, col=3)
typ.names <- c("Ratio of MAD of Bias-Corrected QR \n to MAD of QR")
abline(h=1, col=1, lwd=2, lty=1)
legend(.025, .7, legend = typ.names, lty = c(2), col=c(3), bty="no", lwd=2)

plot( range(.005,.1), range(0,60), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="B. Weibull(1)-disturbance: Estimating Slopes", lwd=2)
lines( taus, volE[seq(2,8, by=2)]/volE[seq(10,16, by=2)], type="l", lty=2, lwd=2, col=3)
typ.names <- c("Ratio of MAD of Bias-Corrected QR \n to MAD of QR")
abline(h=1, col=1, lwd=2, lty=1)
legend(.025, .7, legend = typ.names, lty = c(2), col=c(3), bty="no", lwd=2)


plot( range(.005,.1), range(0,60), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="B. Weibull(4)-disturbance: Estimating Intercept", lwd=2)
lines( taus, volF[seq(1,8, by=2)]/volF[seq(9,16, by=2)], type="l", lty=2, lwd=2, col=3)
typ.names <- c("Ratio of MAD of Bias-Corrected QR \n to MAD of QR")
abline(h=1, col=1, lwd=2, lty=1)
legend(.025, .7, legend = typ.names, lty = c(2), col=c(3), bty="no", lwd=2)

plot( range(.005,.1), range(0,60), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="B. Weibull(4)-disturbance: Estimating Slopes", lwd=2)
lines( taus, volF[seq(2,8, by=2)]/volF[seq(10,16, by=2)], type="l", lty=2, lwd=2, col=3)
typ.names <- c("Ratio of MAD of Bias-Corrected QR \n to MAD of QR")
abline(h=1, col=1, lwd=2, lty=1)
legend(.025, .7, legend = typ.names, lty = c(2), col=c(3), bty="no", lwd=2)

dev.off();