# Step 0 - grab data, let's call it "x" and "y"

xx<-ss.soi$V1
yy<-ss.recruit$V1

# Do the series require detrending?

kdescripb(x,"SS SOI")  
kdescripb(y,"SS Recruit")

# No detrending with the exception of subtracting the mean

y<-yy-mean(yy)
x<-xx-mean(xx)

# STEP 1 - We can capture the main structure in SOI by an AR(12)

xfit<- arima.mle(x,model=list(order=c(14,0,0)))
xfitdiag<-arima.diag(xfit,resid=T)

# STEP 2 - Transform y and x

ytr<-filter(y,filter=c(1,-1*xfit$model$ar),sides=1)
xtr<-filter(x,filter=c(1,-1*xfit$model$ar),sides=1) #equal to xfitdiag$resid

# STEP 3 - Look at the cross correlation between the two transformed series

par(mfrow=c(1,1))
acf(cbind(xtr,ytr))
title(xlab="Cross Correlation from Step 3, the two transformed series")

# Based on the cross correlation, we decide the appropriate lag is 5
# and w(B)=1-wB

# STEP 4 - Build the appropriate regression model

# y(t)=w*y(t-1)+d*x(t-5)+u(t) 
# Fit this regression using OLS to obtain estimates for w and d

n<-length(y)
X<-cbind(y[2:(n-4)],x[6:n])
bivfit<-lsfit(X,y[1:(n-5)],intercept=F)

# STEP 5 - Find and model the noise term eta

# First find the inverse of the MA coefficients
# 1 + wB + (wB)^2 + ...

wma<-bivfit$coef[1]
for (i in 2:50) wma<-c(wma,wma[1]^i)

# Now apply filter to the residual process

eta<-filter(bivfit$resid,c(1,wma),sides=1)

kdescripb(eta,"Additive Noise Process")

# Examining the pacf an AR(2) appears to be the appropriate model
# for eta

etafit<-ar(eta,order.max=2,aic=F,method="burg")
kdescripb(etafit$resid,"Residuals from AR(2) Fit to Additive Noise")  
# you might worry a bit about this fit???


# STEP 6 - Check model fits

# Look at residuals from STEP 1 - we did this
# Look at the residuals from STEP 5 - we did this

# Look at the cross-correlations between the two residual series

par(mfrow=c(1,1))
acf(cbind(xfit$resid[6:n],etafit$resid))