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

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

par(mfrow=c(2,1))
tsplot(soi)
tsplot(rec)

par(mfrow=c(1,1))

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

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

par(mfrow=c(3,4))

for(h in 1:12){
	plot(ts.intersect(rts(soi),lag(rts(soi),h)))
}

par(mfrow=c(3,4))

for(h in 1:12){
	plot(ts.intersect(rts(rec),lag(rts(rec),h)))
}

par(mfrow=c(3,4))

for(h in 1:12){
	plot(ts.intersect(rts(soi),lag(rts(rec),h)))
}

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

x1 <- rts(soi)
x2 <- rts(rec)
u <- ts.intersect(lag(x2,6), x1)   #col 1 of u is rec(t+6), col 2 is soi(t)
cor(u)
          
fit <- lm(u[,1]~u[,2])
summary(fit)

# or

x1 <- rts(soi)
x2 <- rts(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)

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

x1 <- rts(soi)
x2 <- rts(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)

# smoothing rec

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

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))

par(mfrow=c(3,1)) 
plot(time, rec, main="smoothing splines")
lines(smooth.spline(time,rec),spar=.0000001)
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)

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