## ar_test.R -- to accompany Question 12 of the second ## example sheet for Part III/Mphil Monte Carlo Inference, ## Lent 2008, Mathematical Tripos ## ## by Robert B. Gramacy, Statistical Laboratory ## University of Cambridge ## read in the AR sampling functions source("ar.R") ## reading in the data and plotting it x <- ts(scan("ar3.txt", quiet=TRUE)) lx <- length(x) plot(x, type="l", main="randomly generated AR(3) data") ## Making the Xk matrix (or matrices) XKs <- make.XKs(x, 4) ## getting the MCMC samples theta.3 <- ar.gibbs(1000, x[5:lx], Xk=XKs[[3]]) ## summarizing the output names(theta.3) summary(theta.3) # 95% interval quantile(theta.3[,1],probs=c(0.025,0.5,0.975)) plot(theta.3\$a.1, type="l", main="MCMC samples of a1") ## one--step ahead predictive means x501.mean <- as.matrix(theta.3[,1:3]) %*% x[lx:(lx-2)] summary(x501.mean) ## one--step ahead predictive distribution (i.e., forecast) x501 <- rnorm(length(x501.mean), x501.mean, theta.3\$se2) x501q <- quantile(x501, c(0.025, 0.5, 0.975)) x501q ## plotting the forecast last <- (lx-40):lx plot(x=last, y=x[last], type="l", xlim=c(lx-40, lx+1), range(x), main="Posterior predictive interval for x501") points(x=501, y=x501q[2], col=2) arrows(x0=c(501,501),y0=c(x501q[c(2,2)]), x1=c(501,501), y1=x501q[c(1,3)], angle=90, col=2) ## fitting AR(2) and AR(4) models theta.2 <- ar.gibbs(1000, x[5:lx], Xk=XKs[[2]]) summary(theta.2) theta.4 <- ar.gibbs(1000, x[5:lx], Xk=XKs[[4]]) summary(theta.4) ## significance of AR(4) model? quantile(theta.4\$a.4, 0.975)