#### try different models and plot RSS vs p+q

rss <- rep(0,9)  ## creates a vector of 20 zeros.
aic <- rep(0,9)
y <- arima0(x2,order=c(0,0,0))  ## first start with this model.
rss[1] <- sum((y$resid)^2)
aic[1] <- y$aic
for(i in c(1:4)){
y <- arima0(x2,order=c((i-1),0,i))
rss[2*i] <- sum((y$resid)^2)
aic[2*i] <- y$aic
y <- arima0(x2,order=c(i,0,i))
rss[2*i+1] <- sum((y$resid)^2)
aic[2*i+1] <- y$aic
cat(i)  ## prints i on the screen, so you see how long it's taking
}
#postscript("15to16.ps")
par(mfrow=c(1,2))
par(omi=c(1.3,.1,.1,.5))
par(pin=c(2,3))
plot(c(0:8),rss,type="l",xlab="p+q",ylab="Residual sum of squares")
mtext(s=3,l=2,"Figure 15")
plot(c(0:8),aic,type="l",xlab="p+q",ylab="AIC")
mtext(s=3,l=2,"Figure 16")
#dev.off()
## here, AIC suggests ARIMA(0,0,1), but RSS suggests ARIMA(2,0,2).

y <- arima0(x2,order=c(0,0,1))
y2 <- y$resid
beta1 <- y$coef[1]  ### beta_1-hat = .39229814



### Original time-plot & prediction errors with ARIMA(0,1,1) model fit.
l3 <- length(x)
afitres <- rep(0, l3)
afitres[1] <- 0
for(i in c(2:l3)){
afitres[i] <- x[i] - x[i-1] - beta1*afitres[i-1]
}
afit <- x - afitres
sqrt(var(afitres))  #RMS error = .1885 
#postscript("17.ps")
par(mfrow=c(1,3))
par(omi=c(1.3,.1,.1,.5))
par(pin=c(2,3))
plot(c(0,l3)/12,range(x),type="n",xlab="",ylab="")
lines(c(1:l3)/12,x)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"yield (data)")
plot(c(0,l3)/12,range(x),type="n",xlab="",ylab="")
lines(c(1:l3)/12,afit)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"yield [ARIMA(0,1,1) fit]")
mtext(side=3,line=2,"Figure 17: Data, ARIMA(0,1,1) fit, & prediction errors")
plot(c(0,l3)/12,range(afitres),type="n",xlab="",ylab="")
lines(c(1:l3)/12,afitres)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"prediction errors (%)")
#dev.off()


### Original time-plot with ARIMA(0,1,0) model.


bfitres <- rep(0, l3)
bfitres[1] <- 0
bfitres[2:l3] <- x[2:l3]-x[1:(l3-1)]
bfit <- x - bfitres
sqrt(var(bfitres))  #RMS error = .2030
#postscript("bfit1.ps")
par(mfrow=c(1,3))
par(omi=c(1.3,.1,.1,.5))
par(pin=c(2,3))
plot(c(0,l3)/12,range(x),type="n",xlab="",ylab="")
lines(c(1:l3)/12,x)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"yield (data)")
plot(c(0,l3)/12,range(x),type="n",xlab="",ylab="")
lines(c(1:l3)/12,bfit)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"yield [ARIMA(0,1,0) fit]")
mtext(side=3,line=2,"Figure 18: Data, ARIMA(0,1,0) fit, & prediction errors")
plot(c(0,l3)/12,range(bfitres),type="n",xlab="",ylab="")
lines(c(1:l3)/12,bfitres)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"prediction errors (%)")
#dev.off()


### Show both prediction errors and fits
l4 <- l1+1
#postscript("2fits.ps")
par(mfrow=c(1,2))
par(omi=c(1.3,.1,.1,.5))
par(pin=c(2,3))
plot(c(l4,l3)/12,range(c(afitres[l4:l3],bfitres[l4:l3],cfitres[l4:l3])),
type="n",xlab="",ylab="")
lines(c(l4:l3)/12,afitres[l4:l3],lty=1)
lines(c(l4:l3)/12,bfitres[l4:l3],lty=2)
mtext(s=1,l=2,"years")
mtext(s=2,l=2,"Prediction errors (% yield)")
mtext(side=3,line=2,
"Figure 19: Prediction errors, solid = ARIMA(0,1,1).  dotted = ARIMA(0,1,0).")
plot(x,afit,xlab="actual values",ylab="predicted values",type="n")
lines(x,afit,lty=2)
lines(x,bfit,lty=3)
abline(a=0,b=1)
mtext(side=3,line=2,
"Figure 20. dashed = ARIMA(0,1,1). dotted = ARIMA(0,1,0).")
#dev.off()

### RMS prediction errors:
sqrt(var(afitres[l4:l3]))  # ARIMA(0,1,1) error = .184
sqrt(var(bfitres[l4:l3]))  # ARIMA(0,1,0) error = .206

### RMS Errors on 1st 17 years of data:
sqrt(var(afitres[1:l4]))  # ARIMA(0,1,1) error = .193
sqrt(var(bfitres[1:l4]))  # ARIMA(0,1,0) error = .205