source("http://www.stat.pitt.edu/stoffer/tsa2/Rcode/itall.R") soi0 = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/soi.dat")) soi0 = ts(scan("soi.dat")) rec0 = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/recruit.dat")) rec0 = ts(scan("recruit.dat")) # Detrend both series: soi.detrended = lm(soi0~time(soi0))$resid soi = ts(soi.detrended) rec.detrended = lm(rec0~time(rec0))$resid rec = ts(rec.detrended) soi.acf = acf2(soi, maxlag=30) win.graph() soi.arima = sarima(soi,p=1, d=0, q=0); soi.arima w.tilde = soi.arima$fit$residuals win.graph() par(mfrow=c(2,1)) acf(w.tilde, lag.max=30) pacf(w.tilde, lag.max=30) win.graph() phi = soi.arima$fit$coef[1] print(phi) y.tilde = filter(rec, c(1, -phi), sides=1) y.tilde = y.tilde[-1] ccf(y.tilde, w.tilde, lag.max=30) library(dynlm) out = dynlm(rec~L(rec,1) + L(soi, 5) - 1) print(summary(out)) omega = out$coef[1] #.8479 u = out$resid eta = filter(u, filter = omega, "recursive") win.graph() par(mfrow=c(2,1)) acf(eta, lag.max=30) pacf(eta, lag.max=30) # Suggests AR(2) win.graph() eta.arima = sarima(eta,p=2, d=0, q=0); eta.arima z = eta.arima$fit$residuals win.graph() z.acf = acf2(z)