##exmaple 3.35 and 3.36

gnp96<-read.table(file="W:/teaching/stat481581/classnotes/chap3/gnp96.dat")
gnp<-ts(gnp96[,2],start=1947, frequency=4)

par(mfrow=c(2,2))
plot(gnp)
acf(gnp,50)
gnpdiff<-diff(gnp) ##first difference of gnp
plot.ts(gnpdiff)
gnpgr<-diff(log(gnp)) ##growth rate
plot.ts(gnpgr)


par(mfrow=c(2,2))
acf(gnpgr,24)
pacf(gnpgr,24)
gnp2diff<-diff(gnpdiff)
plot.ts(gnp2diff)


##Arima fits
gnpgr.ar<-arima(gnpgr,order=c(1,0,0))
gnpgr.ar
gnpgr.ma<-arima(gnpgr,order=c(0,0,2))
gnpgr.ma
ARMAtoMA(ar=.35,ma=0,10)##prints psi-weights


tsdiag(gnpgr.ma, gof.lag=20) ##diagnostics
hist(gnpgr.ma$resid, br=12)  ##breaks: the n+1 cell boundaries 
qqnorm(gnpgr.ma$resid)
qqline(gnpgr.ma$resid)
shapiro.test(gnpgr.ma$resid)


tsdiag(gnpgr.ar, gof.lag=20) ##diagnostics
par(mfrow=c(2,1))
hist(gnpgr.ar$resid, br=12)  ##breaks: the n+1 cell boundaries 
qqnorm(gnpgr.ar$resid)
qqline(gnpgr.ar$resid)
shapiro.test(gnpgr.ar$resid)

##choose the final model
##AIC
gnpgr.ma$aic
gnpgr.ar$aic

##AICc
log(gnpgr.ma$sigma2)+(222+2)/(222-2-2)
log(gnpgr.ar$sigma2)+(222+1)/(222-1-2)


##sic OR BIC
log(gnpgr.ma$sigma2)+(2*log(222)/222)
log(gnpgr.ar$sigma2)+(1*log(222)/222)


##example 3.37
varve<-read.table(file="C:/jenn/teachingrelated/stat481581/chap3/varve.dat"))
#varve<-read.table(file="W:/teaching/stat481581/classnotes/chap3/varve.dat"))
varve.ma<-arima(log(varve), order=c(0,1,1))
varve.ma ##to display results
tsdiag(varve.ma)

varve.arma<-arima(log(varve), order=c(1,1,1))
varve.arma ##to display results
tsdiag(varve.arma, golf.lag=20)

##choose the final model
##AIC
varve.ma$aic
varve.arma$aic

##AICc
log(varve.ma$sigma2)+(222+2)/(222-2-2)
log(varve.arma$sigma2)+(222+1)/(222-1-2)


##sic OR BIC
log(varve.ma$sigma2)+(2*log(222)/222)
log(varve.arma$sigma2)+(1*log(222)/222)


##example 3.41

phi<-c(rep(0,11),.8)
acf<-ARMAacf(ar=phi,ma=-.5,50)   ## ARMAacf Compute the theoretical autocorrelation function or partial
                                 ## autocorrelation function for an ARMA process.

pacf<-ARMAacf(ar=phi,ma=-.5,50,pacf=T)
par(mfrow=c(1,2))
plot(acf, type="h",xlab="lag")
abline(h=0)
plot(pacf,type="h",xlab="lag")
abline(h=0)


##example 3.43

prod<-read.table(file="C:/jenn/teachingrelated/stat481581/chap3/prod.dat")
par(mfrow=c(3,1))
plot.ts(prod)
acf(prod, 48)
pacf(prod,48)

par(mfrow=c(3,1)) ##p(acf) of d1 data
plot.ts(diff(prod$V1))
acf(diff(prod$V1), 48)
pacf(diff(prod$V1),48)

par(mfrow=c(2,1)) ##p(acf) of d1-d12 data
acf(diff(diff(prod$V1),12), 48)
pacf(diff(diff(prod$V1),12),48)

##fit model
prod.fit3<-arima(prod,order=c(1,1,1),seasonal=list(order=c(2,1,1),period=12))
prod.fit3 ##to view results
tsdiag(prod.fit3, gof.lag=48) #diagnostics

###Forecasts for the final model
prod.pr<-predict(prod.fit3,n.ahead=12)
U<-prod.pr$pred + 2*prod.pr$se
L<-prod.pr$pred - 2*prod.pr$se
month<-337:372
plot(month,prod$V1[month],type="o",xlim=c(337,384),ylim=c(100,180),ylab="production")
lines(prod.pr$pred,col="red",type="o")
lines(U,col="blue",lty="dashed")
lines(L,col="blue",lty="dashed")
abline(v=372.5,lty="dotted")