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


## Here the impulse-response analysis of SOI and Recruits is carried out in two ways

help(LagReg)

########### First method:  Treat SOI as input and Recruits as output.

# The function LagReg does impuse-response analysis of two series and plots the results.
# Output includes the M-1 symmetric lags 's' and the coefficients beta(s) for |s| < M/2.

# To use it, run "out = impulse.response(input, output, L, M, threshold, inverse)". See the help file.

ccf(rec,soi)

win.graph()
out.ir.1 = LagReg(input = soi, output = rec, L=15, M=32, threshold=6.9)

########### Second method:  Treat SOI as output and Recruits as input, then re-arrange.

win.graph()
out.ir.2 = LagReg(rec, soi, L=15, M=32, inverse=TRUE, threshold=.005)
	# NOTE: Get the lags from the plot, not the printed output

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

library(dynlm)  
fit3 = dynlm(rec~L(rec,1) + L(soi,5))
fit3


pred2 = 19.2483 + .7515*lag(rec,1) -47.1695*lag(soi,5)
pred3 = 11.3136+.8434*lag(rec,1) -20.3004*lag(soi,5)
win.graph()
par(mfrow=c(2,1))
ts.plot(rec,pred2,lty=2:1, col=1:2, 
  main = "REC (broken line) and predictions")
MSE2 = sum((rec-pred2)^2)/453 
ts.plot(rec,fit3$fitted,lty=2:1, col=1:2, 
  main = "REC (broken line) and predictions") 
MSE3 = sum((rec-fit3$fitted)^2)/453