# The first time you use R (and perhaps after that as well) you will want to enter these three commands:
#install.packages("astsa")
require(astsa)
help(package = astsa) 

spots = sunspotz # The astsa sereis; 'sunpots' is the inbuilt R series
win.graph()
plot(spots)

## Compute the sample acf
win.graph()
qwe=acf(spots ,36)   
mat = cbind(0:36, qwe$acf)  # 'Bind' these two columns into a matrix
colnames(mat) = c("lag", "acf")
mat
var(spots )

## Compute the p-value for the lag-6 estimated autocorrelation of the "spots" series:
p = 2*(1-pnorm(.5248))  # pnorm(z) gives the probability that Z < z, when Z is N(0,1)
p

## See what differencing does:
spots[1:5]
diff(spots)[1:4]
diff(spots, lag = 2)[1:3]
diff(spots, diff = 2)[1:3]

## Analyze the SOI and Recruits data
ma = filter(soi, sides=2, c(.5, rep(1,11), .5)/12)
win.graph()
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")
	# ts.plot plots several time series on the same axes
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 ~ time(x) + 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)