#  Code for simulating the autocorrelation from the four models
# (slide  18 from The Statistical Model class)

library(MASS)

y <- (nhtemp - 32) / 1.8
nsim <- 10000
ris <- matrix(NA,nrow=nsim, ncol=4)
m <- mean(y)
s <- sd(y)

mle<- fitdistr(y,"t",df=5)
mt <- mle$est[1]
st <- mle$est[2]
x <- 1912:1971-1
mod <- lm(y~x)
smod <- sqrt(sum(resid(mod)^2)/mod$df)
fit2 <- arima(y,xreg=x,order=c(1,0,0))
mu2 <- fit2$coef[2]+fit2$coef[3]*x
for(i in 1:nsim)
{
  ysim <- rnorm(60, m, s)
  ris[i,1] <- as.numeric(unlist(acf(x=ysim, plot=FALSE)[1])[1])
  ysim <- rt(60, df=5)*st+mt
  ris[i,2] <- as.numeric(unlist(acf(x=ysim, plot=FALSE)[1])[1])
  ysim <- rnorm(60, 0, smod) + fitted(mod)
  ris[i,3] <- as.numeric(unlist(acf(x=ysim, plot=FALSE)[1])[1])
  e <- arima.sim(model=list(ar=fit2$coef[1]),n=60, sd=sqrt(fit2$sigma2))
  ysim <- mu2 + e
  ris[i,4] <- as.numeric(unlist(acf(x=ysim, plot=FALSE)[1])[1])
}
boxplot(ris[,1:4], names=c("Model 1","Model 2","Model 3","Model 4"))
abline(h=as.numeric(unlist(acf(x=y, plot=FALSE)[1])[1]), col=2)