# fit recruits

rec_scan("E:\\Stat207\\Class\\Examples\\Ch2\\recruit.dat")
x <- rts(rec)   # x is the time series object
par(mfrow=c(3,1))
ts.plot(x, main="Recruits (x)")
acf(x)
acf(x,type="partial")

#indicates AR(2) fit  but we'll also fit an ARMA(1,1)
# first Yule-Walker to fit AR(2)
rec.ar2 <- ar(x, order.max = 2)
# type 
help(ar) #to see what you get, e.g.
rec.ar2$ar  #gives parameters
rec.ar2$aic # gives AIC from lag 0 to order.max
# now the MLE 
help(arima.mle) # for details
# NOTE arima.mle DOES NOT FIT A CONSTANT - IT ASSUMES MEAN=0
x <- x-mean(x)  #subtract mean
rec.ar2mle <- arima.mle(x, model=list(order=c(2,0,0)), n.cond = 2)
# some results
rec.ar2mle$model
rec.ar2mle$var.coef
rec.ar2mle$aic
rec.ar2mle$sigma2

# Diagnostics
par(mfrow=c(1,1))
arima.diag(rec.ar2mle)

# now the ARMA(1,1) fit
rec.arma <- arima.mle(x, model=list(order=c(1,0,1)), n.cond = 2)
rec.arma$model
rec.arma$var.coef
rec.arma$aic

# Forecast 12 time units ahead for AR(2) model
rec.ar2.fore <- arima.forecast(x,n=12,model=rec.ar2mle$model)
cbind(rec.ar2.fore$mean,rec.ar2.fore$std.err)
par(mfrow=c(1,1))
ts.plot(x,rec.ar2.fore$mean,rec.ar2.fore$mean+2*rec.ar2.fore$std.err,
	rec.ar2.fore$mean-2*rec.ar2.fore$std.err)