### This is S-Plus code for the SOI data. Splus

library(ts)

library(modreg)


soi <- scan("E:\\Stat207\\Class\\Examples\\Ch1\\soi.dat")
rec <- scan("E:\\Stat207\\Class\\Examples\\Ch1\\recruit.dat")

soi <- scan("http://www.sci.csueastbay.edu/~esuess/Statistics_6871/Handouts/Examples/Ch1/soi.dat")
rec <- scan("http://www.sci.csueastbay.edu/~esuess/Statistics_6871/Handouts/Examples/Ch1/recruit.dat")

win.graph()

par(mfrow=c(2,1))
ts.plot(soi)
ts.plot(rec)

win.graph()

acf(cbind(soi,rec))  # ACFs and CCF

win.graph()

pairs.default(cbind(soi,rec))  # scatter plot matrix

# Regress one series on another [rec(t) on soi(t-6)]

x1 <- ts(soi)
x2 <- ts(rec)
u <- ts.intersect(lag(x2,6), x1)   #col 1 of u is rec(t+6), col 2 is soi(t)

plot(u[,1],u[,2])
cor(u)
          
fit <- lm(u[,1]~u[,2])
summary(fit)

noise <- fit$residuals  # get residuals  

win.graph()

acf(noise)

win.graph()

ts.plot(noise)

### This is S-Plus code for the SOI data. R

win.graph()

par(mfrow=c(2,1))
ts.plot(soi)
ts.plot(rec)

win.graph()

par(mfrow=c(1,1))

acf(cbind(soi,rec))  # ACFs and CCF

win.graph()

lag.plot(soi, lag=12, layout=c(3,4))  # lag plots

win.graph()

lag.plot(rec, lag=12, layout=c(3,4))  

# Regress one series on another [rec(t) on soi(t-6)]

x1 <- ts(soi)
x2 <- ts(rec)
u <- ts.intersect(x2, lag(x1,-6))   #col 1 of u is rec(t), col 2 is soi(t-6)

cor(u)
          
fit <- lm(u[,1]~u[,2])
summary(fit)

noise <- fit$residuals  # get residuals  

win.graph()

acf(noise)

win.graph()

ts.plot(noise)

# Regress one series on the past of the the other [rec(t) on soi(t-1),...,soi(t-8)]

x1 <- ts(soi)
x2 <- ts(rec)
u <- ts.intersect(x2, lag(x1,0), lag(x1,-1), lag(x1,-2), lag(x1,-3), lag(x1,-4), lag(x1,-5), lag(x1,-6), lag(x1,-7), lag(x1,-8))   #col 1 of u is rec(t+6), col 2 is soi(t)

cor(u)

fit <- lm(u[,1]~u[,2:10])
summary(fit)

noise <- fit$residuals  # get residuals  

win.graph()

acf(noise)

win.graph()

ts.plot(noise)

# smoothing rec

time <- c(1:length(rec))

win.graph()

par(mfrow=c(3,1)) 
plot(time, rec, main="moving average")
filter5 <-filter(rec, c(1,1,1,1,1)/5)
lines(time, filter5)
filter53 <- filter(rec, rep(1,53)/53 )
lines(time, filter53)

plot(time, rec, main="nearest neighbor")
lines(supsmu(time, rec,bass=10))
lines(supsmu(time, rec,span=.01))

plot(time, rec, main="ksmooth")
lines(ksmooth(time,rec,"normal",bandwidth=4))
lines(ksmooth(time,rec,"normal",bandwidth=100))

win.graph()

par(mfrow=c(3,1)) 
plot(time, rec, main="smoothing splines")
lines(smooth.spline(time,rec))
lines(smooth.spline(time,rec,spar=.1))

plot(time, rec, main="natural splines")
lines(smooth.spline(time,rec))
lines(smooth.spline(time,rec,spar=0,df=10))

plot(time, rec, main="lowess")
lines(lowess(time,rec,.02))
lines(lowess(time,rec))

# rec versus lag(soi,-6)

win.graph()

plot(lag(soi,-6), rec)  
lines(lowess(lag(soi,-6), rec,.02))
lines(lowess(lag(soi,-6), rec))