flu = scan("http://www.stat.ualberta.ca/~wiens/stat479/S&Sdatasets/flu.dat")
flu = as.ts(flu)
par(mfrow=c(3,1))
plot(flu, type="o",ylab="deaths per 10,000", xlab="", xaxp=c(12,132,10))
plot(rev(flu), type="o",ylab="deaths in reverse order", 
	xlab="", xaxp=c(12,132,10))

x0 = diff(flu)

library(dynlm)  
# This function dynlm fits linear models with lagged series; 
# it addresses the problem that the lagged series will have differing lengths
out0 = dynlm(x0 ~ L(x0,1) + L(x0,2) + L(x0,3) + L(x0,4), x=T, y=T) 
 # Fit this model (with x=T, y=Y) merely to get the right columns of lagged variables
model = out0$model # These will be here
model
x = model[,1]
X = model[,2:5]
mean(x[x>0]) ## = .052
abs(mean(x[x<0]))  ## =.058; flu decreases more quickly than it increases
plot(x, type="o",ylab="diff(flu)", xlab="", xaxp=c(12,132,10))
abline(h=.05)

less = (X[,1]<.05)
x1 = x[less]
X1 = X[less,]
colnames(X1) = c("x(t-1)","x(t-2)","x(t-3)","x(t-4)")
out1 = lm(x1~X1[,1]+X1[,2]+X1[,3]+X1[,4]) ##S&S fit intercept, but don't report it
print(summary(out1))

greater = (X[,1]>=.05)
x2 = x[greater]
X2 = X[greater,]
colnames(X2) = c("x(t-1)","x(t-2)","x(t-3)","x(t-4)")
out2 = lm(x2~X2[,1]+X2[,2]+X2[,3]+X2[,4]) ##Here they do report it
print(summary(out2))


qwe1=out1$fitted.values
qwe2 = predict(out1)
cbind(qwe1,qwe2)


fit1 = predict(out1)
fit2 = predict(out2)
less[less==1]= fit1
greater[greater==1] = fit2
fit = less + greater
win.graph()
plot(x, type="o",ylab="diff(flu)", xlab="", xaxp=c(12,132,10))
lines(fit, col = "blue", lty="dashed")