#The MC is done in S-Plus. Graphics is done in R.
#Y is generated according to a t-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 ####################################

inqr<- function(Y) { return(quantile(Y, .1), quantile(Y, .9))}


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



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



filenA<-  "c:\\aaa\\comput\\CEC\\nn500.CASEA.HS.run1.txt"
filen.biasA<-  "c:\\aaa\\comput\\CEC\\nn500.CASEA.HS.BiasAnalysis.run1.txt"
filenB<-  "c:\\aaa\\comput\\CEC\\nn500.CASEB.HS.run1.txt"
filen.biasB<-  "c:\\aaa\\comput\\CEC\\nn500.CASEB.HS.BiasAnalysis.run1.txt"
filenC<-  "c:\\aaa\\comput\\CEC\\nn500.CASEC.HS.run1.txt"
filen.biasC<-  "c:\\aaa\\comput\\CEC\\nn500.CASEC.HS.BiasAnalysis.run1.txt"



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

datt<- gen.data(nn, DF=3);  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/5}  # A Hall-Sheather Rule


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

# generate data ; selct badnwidth

datt<- gen.data(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*qt(.005,DF),  extr.005$BC.e[3] - 1);
mise.020<- c(  extr.020$BC.e[1] - 1- .07*qt(.020,DF),  extr.020$BC.e[3] - 1);
mise.050<- c(  extr.050$BC.e[1] - 1- .07*qt(.010,DF),  extr.050$BC.e[3] - 1);
mise.100<- c(  extr.100$BC.e[1] - 1- .07*qt(.100,DF),  extr.100$BC.e[3] - 1);
 

cove.005<-  c( I( extr.005$ciL.e[1] < 1+ .07*qt(.005,DF) )* I ( 1+ .07*qt(.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*qt(.020,DF) )* I ( 1+ .07*qt(.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*qt(.050,DF) )* I ( 1+ .07*qt(.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*qt(.100,DF) )* I ( 1+ .07*qt(.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*qt(.005,DF) )* I ( 1+ .07*qt(.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*qt(.020,DF) )* I ( 1+ .07*qt(.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*qt(.050,DF) )* I ( 1+ .07*qt(.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*qt(.100,DF) )* I ( 1+ .07*qt(.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*qt(.005,DF) ,  centr.005$betat[3] - 1 );
misc.020<- c(  centr.020$betat[1] - 1- .07*qt(.020,DF) ,  centr.020$betat[3] - 1 );
misc.050<- c(  centr.050$betat[1] - 1- .07*qt(.010,DF) ,  centr.050$betat[3] - 1 );
misc.100<- c(  centr.100$betat[1] - 1- .07*qt(.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:250) {monte.cov(filenA, filen.biasA, DF=1)} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:250) {monte.cov(filenB, filen.biasB, DF=3)} ; speed <- (proc.time()-now)/3600;
now<- proc.time(); for(i in 1:250) {monte.cov(filenC, filen.biasC, DF=30)} ; speed <- (proc.time()-now)/3600;


# This part is DONE in R (!!!):

volA<- apply(read.table(filenA),2, mean)
volB<- apply(read.table(filenB),2, mean)
volC<- apply(read.table(filenC),2, mean)

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


postscript("c:\\aaa\\\comput\\cec\\MCresultABC.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="A. t(1)-disturbance: Coverage for Intercept", lwd=2)
lines( taus, volA[seq(1,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volA[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="A. t(1)-disturbance: Coverage for Slopes", lwd=2)
lines( taus, volA[seq(2,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volA[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="B. t(3)-disturbance: Coverage for Intercept", lwd=2)
lines( taus, volB[seq(1,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volB[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="B. t(3)-disturbance: Coverage for Slopes", lwd=2)
lines( taus, volB[seq(2,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volB[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="C. t(30)-disturbance: Coverage for Intercept", lwd=2)
lines( taus, volC[seq(1,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volC[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="C. t(30)-disturbance: Coverage for Slopes", lwd=2)
lines( taus, volC[seq(2,8, by=2)], type="l", lty=1, lwd=2, col=4)
lines( taus, volC[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();






plot(apply(read.table(filenA),2, mean))


apply(read.table(filenB),2, mean)
apply(read.table(filenC),2, mean)


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

volA<- apply(read.table(filen.biasA),2, medad)
volB<- apply(read.table(filen.biasB),2, medad)
volC<- apply(read.table(filen.biasC),2, medad)




postscript("c:\\aaa\\\comput\\cec\\MC.Forecast.resultABC.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.2)

plot( range(.005,.1), range(.5,1.5), type="n", xlim=xlim, ylim=ylim, ylab=ylab, xlab=xlab, main="A. t(1)-disturbance: Estimating Intercept", lwd=2)
lines( taus, volA[seq(1,8, by=2)]/volA[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. t(1)-disturbance: Estimating Slopes", lwd=2)
lines( taus, volA[seq(2,8, by=2)]/volA[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. t(3)-disturbance: Estimating Intercept", lwd=2)
lines( taus, volB[seq(1,8, by=2)]/volB[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. t(3)-disturbance: Estimating Slopes", lwd=2)
lines( taus, volB[seq(2,8, by=2)]/volB[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. t(30)-disturbance: Estimating Intercept", lwd=2)
lines( taus, volC[seq(1,8, by=2)]/volC[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. t(30)-disturbance: Estimating Slopes", lwd=2)
lines( taus, volC[seq(2,8, by=2)]/volC[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();






tsplot(abs(mmmB[1:8])/ abs(mmmB[9:16]))