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

prod_rts(scan("E:\\Stat207\\Class\\Examples\\Ch2\\prod.dat"))
ts.plot(prod)		# nonstationary
acf(prod)
acf(prod,type="partial")

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

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

# Try p = 1

model1 <- list(list(order=c(1,1,0)), list(order=c(0,1,0), period=12))
fit1_arima.mle(prod, model=model1)
fit1
arima.diag(fit1)		# 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_arima.mle(prod, model=model2)
fit2
arima.diag(fit2)