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