################################################################# ### Example 7.11: Bootstrap intervals for quantiles ################################################################# ### # Functions ### quantestsample<-function(size,samplesize,p,S0,mu,sigma) { vals<-(1:size) for(i in (1:size)) { vals[i]<-empquant(S0*(1-rlnorm(samplesize,mu,sigma)),0.95) } vals } empquant<-function(datavec,p) { sort(datavec,decreasing=TRUE)[floor(length(datavec)*(1-p))+1] } empquantboot<-function(datavec,indvec,p) sort(datavec[indvec],decreasing=TRUE)[floor(length(datavec)*(1-p))+1] mybootdistr<-function(datavec,rep,p) { vals<-(1:rep) datasize<-length(datavec) for(i in (1:rep)) { vals[i]<-empquantboot(datavec,sample((1:datasize),datasize,replace=TRUE),p) } vals } mybootci<-function(datavec,estim,rep,p,q) { vals<-(1:rep) mbci<-(1:2) vals<-estim - mybootdistr(datavec,rep,p) mbci<-c(estim + empquant(vals,(1-q)/2),estim + empquant(vals,(1+q)/2)) mbci } ### # Calls to functions and plots ### nrows<-500 ncols<-500 S0<-1 alpha<-0.05 mu<-0 sigma<-0.01 returnmat<-matrix(0,nrows,ncols) negxvalmat<-returnmat for(i in (1:nrows)) { returnmat[i,]<-rlnorm(ncols,mu,sigma) negxvalmat[i,]<-S0*(1-returnmat[i,]) } quantestvals<-quantestsample(2000,500,0.95,1,0,0.01) hist(0.01631400-quantestvals,100,xlab="",ylab="",main="",yaxt="n",col="black") # centeredquantestsample_2000.eps hist(empquant(negxvalmat[1,],0.95)-mybootdistr(negxvalmat[1,],2000,0.95),100,xlab="",ylab="",main="",yaxt="n",col="black") # centeredbootdistr1_2000_100.eps hist(empquant(negxvalmat[2,],0.95)-mybootdistr(negxvalmat[2,],2000,0.95),100,xlab="",ylab="",main="",yaxt="n",col="black") # centeredbootdistr2_2000_100.eps hist(empquant(negxvalmat[3,],0.95)-mybootdistr(negxvalmat[3,],2000,0.95),100,xlab="",ylab="",main="",yaxt="n",col="black") # centeredbootdistr3_2000_100.eps hist(empquant(negxvalmat[1,],0.95)-mybootdistr(negxvalmat[1,],10000,0.95),100,xlab="",ylab="",main="",yaxt="n",col="black") # centeredbootdistr1_10000_100.eps hist(empquant(negxvalmat[2,],0.95)-mybootdistr(negxvalmat[2,],10000,0.95),100,xlab="",ylab="",main="",yaxt="n",col="black") # centeredbootdistr2_10000_100.eps bootcis<-matrix(0,100,2) for(i in (1:100)) { bootcis[i,]<-mybootci(negxvalmat[i,],empquant(negxvalmat[i,],0.95), 2000,0.95,0.95) } exactcis<-matrix(0,100,2) for(i in (1:100)){ exactcis[i,]<-c(sort(negxvalmat[i,],decreasing=TRUE)[35],sort(negxvalmat[i,],decreasing=TRUE)[16]) } plot((1:50),-(1:50),xlab="",ylab="",ylim=c(min(min(exactcis[1:50,1]),min(bootcis[1:50,1])),max(max(exactcis[1:50,2]),max(bootcis[1:50,2])))) for(i in (1:50)){ segments(i,exactcis[i,1],i,exactcis[i,2]) segments(1,0.01631400,50,0.01631400) } # confintsexactplustrue.eps plot((1:50),-(1:50),xlab="",ylab="",ylim=c(min(min(exactcis[1:50,1]),min(bootcis[1:50,1])),max(max(exactcis[1:50,2]),max(bootcis[1:50,2])))) for(i in (1:50)) segments(i,bootcis[i,1],i,bootcis[i,2]) segments(1,0.01631400,100,0.01631400) # confintsbootplustrue.eps