### Estimate arima models

load("/Users/bob/courses/stat910/rcode/tsa3.rda")   

# --- Simulate ARMA(2,2) ... generates a time series object

ar.coef <- c(1.3,-0.5);
polyroot(c(1,-ar.coef))

ma.coef <- c(0.5,0.2)
polyroot(c(1,ma.coef))


n <- 400;
y <- arima.sim(n=n, model=list(ar=ar.coef, ma=ma.coef), rand.gen=rnorm)

plot(y)


# --- estimation... lots of options, methods
?arima
arma.fit <- arima(y, order=c(2,0,2)); arma.fit


# --- large AR approximation
fit.ar <- arima(y,order=c(20,0,0))               # slow....

?ar
fit.ar <- ar.burg(y, aic=FALSE, order.max=20,)   # fast....
summary(fit.ar)

se <- sqrt(diag(fit.ar$asy.var.coef))
cbind(fit.ar$ar, se, round(abs(fit.ar$ar)/se,2) )

w <- as.ts(fit.ar$resid); plot(w)               # convert into a time series

join <- ts.intersect(y, lagy.1=lag(y,1), lagy.2=lag(y,2), lagw.1=lag(w,1), lagw.2=lag(w,2))

regr <- lm(y ~ lagy.1 + lagy.2 + lagw.1 + lagw.2, data=join)
summary(regr)

join[110:113,]    # Eek!  See S&S rant


# --- fit via regression  (caution on 'intercept')

join <- ts.intersect(y, lagy.1=lag(y,-1), lagy.2=lag(y,-2), lagw.1=lag(w,-1), lagw.2=lag(w,-2))
join[110:113,]   

regr <- lm(y ~ lagy.1 + lagy.2 + lagw.1 + lagw.2, data=join)
summary(regr)

arma.fit