## Analyze the SOI and Recruits data
## Either of these ways of getting the data will work:
# This gets the data from the website:
soi = ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/soi.dat"), start=1950, frequency=12)
rec =  ts(scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/recruit.dat"), start=1950, frequency=12)

# This will get the file from your pc, provided that R has been started from the directory, on your pc, containing the datasets:
soi = ts(scan("soi.dat"), start=1950, frequency=12)
rec = ts(scan("recruit.dat"), start=1950, frequency=12)

# You can specify the directory from which R starts - right click on the R icon and go to 'Properties'.
# Of course YOU will use only one of these two ways of getting the data; the other might give an error message but this can be ignored.

# 'scan' for univariate series, 'read.table' for multivariate

ma = filter(soi, sides=2, c(.5, rep(1,11), .5)/12)

par(mfrow=c(2,1))
ts.plot(soi, ma, lty=2:1, col=1:2, lwd=1:2, main="Southern Oscillation Index and 12 month moving average")
plot(rec, main = "Recruits")

win.graph()
par(mfrow=c(3,1))
acf(soi, 50)
acf(rec, 50)
ccf(soi, rec, 50)



## For this next step, first click on "Packages" in the menu at the top of the R screen.
## Choose the Canadian mirror, and then when prompted ask to have the 'dynlm' package installed.

library(dynlm)  
# This function dynlm fits linear models with lagged series; 
# it addresses the problem that the lagged series will have differing lengths

## Regress Y = Recruits on X = SOI, at various lags.
x = soi - mean(soi)  # Address possible collinearity
trend = time(x)
fit = dynlm(rec ~ trend + L(x,3) + L(x,4)+ L(x,5) + L(x,6) + L(x,7) + L(x,8)+ L(x,9) + L(x,10)+ L(x,11) + L(x,12) )

# Or merely "fit = dynlm(rec ~ trend + L(x,3:12))"

summary(fit)

win.graph()
par(mfrow=c(3,1))
resids = fit$residuals
fits = fit$fitted.values
plot.ts(resids)
plot(fits, resids, main = "Residuals vs. Fitted values")
abline(h=0)
acf(resids)