# Fit an ARIMA(1,1,0)x(0,1,1)12 to the Production Series

library(ts)

#prod <- ts(scan("E:\\Stat207\\Class\\Examples\\Ch2\\prod.dat"))

prod <- scan("http://www.sci.csueastbay.edu/~esuess/Statistics_6871/Handouts/Handout7/PROD.DAT")

win.graph()
ts.plot(prod)		# nonstationary
win.graph()
acf(prod)
win.graph()
acf(prod,type="partial")

prod.diff <- diff(prod)
win.graph()
ts.plot(prod.diff)		# seasonal nonstationary
win.graph()
acf(prod.diff)
win.graph()
acf(prod.diff,type="partial")

prod.diff.sdiff <- diff(prod.diff,12)
win.graph()
ts.plot(prod.diff.sdiff)	
win.graph()
acf(prod.diff.sdiff)		
win.graph()
acf(prod.diff.sdiff,type="partial")		# note the spike at 1 in the PACF, fit p=1

# Try p = 1

fit1 <- arima0(prod, order=c(1,1,0), seasonal=list(order=c(0,1,0), period=12))
summary(fit1)

win.graph()
par(mfrow=c(3,1))
ts.plot(fit1$resid)
acf(fit1$resid)
acf(fit1$resid,type="partial")	# note the spike at lag 12 in the ACF, fit Q=1

model2 <- list(list(order=c(1,1,0)), list(order=c(1,1,0), period=12))
fit2_arima0(prod, order=c(1,1,0), seasonal=list(order=c(1,1,0), period=12))
summary(fit2)

par(mfrow=c(3,1))
ts.plot(fit2$resid)
acf(fit2$resid)
acf(fit2$resid,type="partial")	# note the spike at lag 12 in the ACF, fit Q=1