1#Analysis of the QAR(1) Melbourne Temperature Example 2require(splines) 3if(require(hdrcde)){ 4 if(interactive()){ 5 oldpar <- par(ask = TRUE) 6 x <- maxtemp[-3650] 7 y <- maxtemp[-1] 8 s <- (x<40) #Delete a few (influential, ridiculously hot) days 9 x <- x[s] 10 y <- y[s] 11 z <- seq(10,36,length=100) 12 13 fit <- rq(y~ bs(x,knots=quantile(x,c(.05,.25,.5,.75,.95))), tau = 1:19/20) 14 par(cex=1,pty="s") 15 xlab <- "yesterday's max temperature" 16 ylab <- "today's max temperature" 17 plot(x,y,pch=".",xlab=xlab,ylab=ylab) 18 matlines(z,predict(fit, newdata = data.frame(x = z)), lty = 1) 19 abline(c(0,1),lty=3) 20 title("Melbourne QAR Model") 21 22 taus <- 1:199/200 23 xs <- c(11,16,21,25,30,35) 24 fit <- rq(y~ bs(x,knots=quantile(x,c(.05,.25,.5,.75,.95))), tau = taus) 25 Qy <- predict(fit,newdata = data.frame(x = xs)) 26 par(mfrow = c(2,3)) 27 for(i in 1:length(xs)){ 28 Qyi <- Qy[i,-1] 29 fhat <- akj(Qyi,Qyi,diff(taus), h = 1)$dens 30 xlab <- "today's max temperature" 31 plot(Qyi,fhat,type="l",xlab=xlab,ylab="density") 32 abline(v=xs[i], col="red") 33 title(paste("Yesterday's Temp", format(round(xs[i])))) 34 } 35 par(oldpar) 36 } 37 else warning("Need hdrcde package to get data for this demo") 38 } 39 40