# Time series fit and model checking for unemployment series unemployment <- read.table ("unemployment.dat", header=TRUE) year <- unemployment\$year y <- unemployment\$unemployed.pct # plot the unemployment rate postscript ("c:/books/multilevel/unemployment1.ps", height=3, width=4) par (mar=c(4,4,2,2)) plot (year, y, type="l", ylab="unemployment", xlab="year", yaxs="i", ylim=c(0, max(y)*1.05), yaxt="n", mgp=c(2,.5,0), cex.axis=1.2, cex.lab=1.2) axis (2, c(0,5,10), paste (c(0,5,10), "%", sep=""), mgp=c(2,.5,0), cex.axis=1.2) dev.off() # fit a 1st-order autogregression n <- length (y) print (n) y.lag <- c (NA, y[1:(n-1)]) lm.lag <- lm (y ~ y.lag) display (lm.lag) # simulate replicated datasets (using beta.hat, sigma.hat) b.hat <- beta.hat (lm.lag) s.hat <- sigma.hat (lm.lag) n.sims <- 1000 y.rep <- array (NA, c(n.sims, n)) for (s in 1:n.sims){ y.rep[s,1] <- y[1] for (t in 2:n){ prediction <- c (1, y.rep[s,t-1]) %*% b.hat y.rep[s,t] <- rnorm (1, prediction, s.hat) } } # simulate replicated datasets (full uncertainty) lm.lag.sim <- sim (lm.lag, n.sims) for (s in 1:n.sims){ y.rep[s,1] <- y[1] for (t in 2:n){ prediction <- c (1, y.rep[s,t-1]) %*% lm.lag.sim\$beta[s,] y.rep[s,t] <- rnorm (1, prediction, lm.lag.sim\$sigma[s]) } } # plot the simulated unemployment rate series postscript ("c:/books/multilevel/unemployment2.ps", height=5, width=10) par (mfrow=c(3,5), mar=c(4,4,2,2)) for (s in 1:15){ plot (year, y.rep[s,], type="l", ylab="unemployment", xlab="year", yaxs="i", ylim=c(0, max(y)*1.05), yaxt="n", mgp=c(2,.5,0), cex.axis=1.5, cex.lab=1.5, main=paste ("simulation", s), cex.main=1.5) axis (2, c(0,5,10), paste (c(0,5,10), "%", sep=""), mgp=c(2,.5,0), cex.axis=1.5) } dev.off() # numerical model check test <- function (y){ n <- length (y) y.lag <- c (NA, y[1:(n-1)]) y.lag2 <- c (NA, NA, y[1:(n-2)]) sum (sign(y-y.lag) != sign(y.lag-y.lag2), na.rm=TRUE) } print (test (y)) test.rep <- rep (NA, n.sims) for (s in 1:n.sims){ test.rep[s] <- test (y.rep[s,]) } print (mean (test.rep > test(y))) print (quantile (test.rep, c(.05,.5,.95)))