#install.packages("astsa")
require(astsa)
#help(package = astsa) 


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

par(mfrow=c(2,1))
# lty = 1: solid, =2: dashed.
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()
acf2(soi, 50)
win.graph()
acf2(rec, 50)
win.graph()
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.
## If you've done this once already you shouldn't have to do it again.

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

## Regress Y = Recruits on X = SOI, at various lags.
fit = dynlm(rec ~ time(x) + 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) )
summary(fit)

win.graph() # Figure 1
par(mfrow=c(2,1))
resids = fit$residuals
fits = fit$fitted.values
plot.ts(resids, main = "Residuals from lagged regression on SOI")
abline(h=0)
ts.plot(rec, fits, lty=2:1, col=1:2, lwd=1:2, main="Recruits and estimated regression response")

win.graph() 
acf2(resids)

##Look for a more parsimonious model
win.graph()  # Figure 2
x = diff(soi, 12)
y = diff(rec, 12)
ts.plot(100*x,y,lty=2:1, col=1:2, lwd=1:2, main = "X = 100*differenced SOI (broken line) and Y = differenced REC")
win.graph()
ccf(x,y, main = "COV(X(t+m),Y(t))")
win.graph()
acf2(x)
win.graph()
acf2(y)

y1 = lag(y, -1)
y2 = lag(y, -2)
x10 = lag(x, -10)
fit2 = dynlm(y ~ y1 + y2 + x10 - 1)  # No intercept is to be fitted
z = ts(fit2$resid)
fit2$coef

win.graph() # Figure 3
par(mfrow=c(2,1))
fits2 = fit2$fitted.values
ts.plot(y, fits2, lty=2:1, col=1:2, lwd=1:2, main = "Y (broken line) and fitted values")
plot(z, main = "Residuals Z(t)")

win.graph() # Figure 4
par(mfrow=c(3,1))
plot(z)
acf2(z, 48)

out = matrix(ncol = 7)
for (p in 0:1) { for (q in 0:1) { for (P in 0:1)  { for (Q in 0:1)  {
fits = sarima(z, p, 0, q, P, 0, Q, 12)
out = rbind(out, c(p, q, P, Q,  fits$AIC, fits$AICc, fits$BIC))
}}}}
out = out[-1,]
colnames(out) = c("p", "q", "P", "Q", "AIC", "AICc","BIC")
out


out = round(out,5)
out[out[,5] == min(out[,5])] # min AIC
out[out[,6] == min(out[,6])] # min AICc
out[out[,7] == min(out[,7])] # min BIC

win.graph() # Figure 5
fit3 = sarima(z, 0,0,1, 1,0,1, 12)

win.graph()
sarima.for(z, n.ahead = 12, 0,0,1, 0,0,1, 12)
abline(v = length(z)-.05)


## Determine the roots of the regression characteristic equation
a = c(1, -1.3086448, 0.4198637)
zeroes = polyroot(a)
zeroes
abs(zeroes)



##################################################

#install.packages("astsa")
require(astsa)
#help(package = astsa) 

# Look at the data (recall you were led through this in Asst. 1):
par(mfrow=c(2,1))
plot(varve, sub = "varve series")
plot(log(varve), sub = "log varve")
win.graph()
acf2(log(varve))

win.graph()
x = diff(log(varve))
plot(x, main = "Diff(log(varve)) series")
win.graph()
acf2(x)
# diff(log(varve)) looks stationary (Interpretation?)
# ARMA(p,1) with p=1,2,3?

sarima(log(varve), p=0,d=1,q=1) # A poor fit; compare with Figure 3.18 in the text

win.graph()
fit1 = sarima(log(varve), 1,1,1)
names(fit1)
names(fit1$fit)

# Looks pretty good; can we do much better?
# Try several models:
out = matrix(ncol = 5)
for (p in 0:4) { for (q in 0:2) {
fits = sarima(log(varve), p, 1, q)
out = rbind(out, c(p, 1, q, fits$AIC, fits$AICc, fits$BIC))
}}
out = out[-1,]
colnames(out) = c("p", "q", "AIC", "AICc","BIC")
out

out = round(out,5)
out[out[,3] == min(out[,3])] # min AIC
out[out[,4] == min(out[,4])] # min AICc
out[out[,5] == min(out[,5])] # min BIC


# So p=1 was a good guess - a triumph of the Principle of Parsimony

fit = fit1

# Forecasts for the final model   
par(mfrow=c(1,1))          
logvarve.pr = sarima.for(log(varve), n.ahead=50, p=1,d=1,q=1)
abline(v=634-.05)

#################################
# How well does the model forecast?
# Hold back 'h' observations:
h = 20
n1 = length(varve) - h
logvarve1 = ts(log(varve[1:n1]))
logvarve2 = ts(log(varve[(n1+1):length(varve)]), start = n1+1)

# Predict logvarve2 with an arima(1,1,1) model fitted from the data in logvarve1:
par(mfrow=c(1,1))
logvarve1.pr = sarima.for(logvarve1, n.ahead = h, 1,1,1)
lines(logvarve2)
abline(v=n1+.5)