d.0=c(4,14,19,27,27,41,25,37,43,47,29,23,9,2)
d.1=c(2, 8,13, 9,14,15,13,14,21,22,11, 4,2,0)
d.2=c(4, 6, 8,13,17,13,14,21,21,21, 8, 3,5,0)

d.0.f = filter(d.0,rep(1/window.w,window.w))
     d.1.f = filter(d.1,rep(1/window.w,window.w))
     d.2.f = filter(d.2,rep(1/window.w,window.w))
     
     hr.12 = (d.1.f+d.2.f)/d.0.f
     hr.1 = d.1.f/(d.0.f/2)
     hr.2 = d.2.f/(d.0.f/2)

par(mfrow=c(1,1),mar = c(4,6,1,2))

     plot( ((1:nb)-0.5)*bin.w,hr.12,type="l",
      ylim=c(0,1.5), xlim=c(0,14),
      xlab="",
      ylab="Hazard Ratio (HR)\n(3-year moving window)",
      yaxt="n")
     
     
     NB=1000 # bootstrap samples
     
     for (B in 1:NB) {
     	
     	d.0.B = rpois(nb,d.0)
     	d.1.B = rpois(nb,d.1)
     	d.2.B = rpois(nb,d.2)
     	
     	d.0.B.f = filter(d.0.B,rep(1/window.w,window.w))
        d.1.B.f = filter(d.1.B,rep(1/window.w,window.w))
        d.2.B.f = filter(d.2.B,rep(1/window.w,window.w))
     	
     	hr.B = (d.1.B.f+d.2.B.f)/d.0.B.f
        lines( ((1:nb)-0.5)*bin.w,hr.B,type="l",
         lwd=0.25, col="grey80")
     	
     }
     lines( ((1:nb)-0.5)*bin.w,hr.12,lwd=4,col="black")
     
#lines( ((1:nb)-0.5)*bin.w,hr.1,lwd=1,col="red")
#lines( ((1:nb)-0.5)*bin.w,hr.2,lwd=2,col="red")
     
     
     abline(h=seq(0.2,1.4,0.2), col="lightblue" )
     abline(h=1, lwd=1, col="blue")
     abline(v=1:13, col="lightblue" )
     
     for(Hr in seq(0.2,1.4,0.2)) text(0,
         Hr,toString(Hr),cex=0.9,adj=c(0.5,0.5))
         
     text(6,0.2,
     "1000 bootstrap HR functions\nshown in grey",
     cex=0.75,col="grey65")