#Read in data: Jan 1980 - Oct 1991
Wine<-scan("http://www.math.mcgill.ca//dstephens//545//Data//wine.txt")
n<-length(Wine)
tvec<-c(1:n)


plot(Wine,pch=19)
lines(tvec,Wine)

Wine.fit<-lm(Wine~tvec)

mt.Wine<-Wine.fit$coef[1]+Wine.fit$coef[2]*tvec

Wine.Yt<-Wine-mt.Wine

plot(Wine.Yt,pch=19)
lines(tvec,Wine.Yt)

#################################################
lWine<-log(Wine)
plot(lWine,pch=19)
lines(tvec,lWine)

lWine.fit<-lm(lWine~tvec)

mt.lWine<-lWine.fit$coef[1]+lWine.fit$coef[2]*tvec

lWine.Yt<-lWine-mt.lWine

plot(lWine.Yt,pch=19)
lines(tvec,lWine.Yt)
#################################################


acf(Wine.Yt,main="ACF of Wine residuals")
acf(lWine.Yt,main="ACF of log(Wine) residuals")

################################################
tsl.lWine<-stl(ts(lWine,frequency=12,start=c(1980,1)),"per")

plot(tsl.lWine)

Yt<-ts(lWine,frequency=12,start=c(1980,1))
mt<-tsl.lWine$time.series[,2]
st<-tsl.lWine$time.series[,1]

plot(Yt)
points(Yt,pch=19)

lines(tvec,mt+st,col='red')

Yt.res<-ts(tsl.lWine$time.series[,3],frequency=12,start=c(1980,1))
plot(Yt.res)
points(Yt.res,pch=19)

acf(Yt.res)


##########################################################################

t.v<-c(1:n)
Time<-t.v
C1<-cos(2*pi*t.v/12)
S1<-sin(2*pi*t.v/12)
C2<-cos(2*pi*t.v/6)
S2<-sin(2*pi*t.v/6)

lWine.data<-data.frame(lWine,Time,S1,C1,S2,C2)

lWine.fit<-lm(lWine~Time+S1+C1+S2+C2,data=lWine.data)
lWine.fitted<-lWine.fit$fitted
lWine.coef<-lWine.fit$coef

plot(Yt,ylab="Monthly wine sales",type='l')
points(Yt,pch=19)
title("Monthly wine sales: log scale \n Jan 1980 - Oct 1991")
lWineFit<-ts(lWine.fitted, frequency = 12, start = c(1980, 1))
points(lWineFit,pch=3,col='red')
lines(lWineFit,lty=2,col='red')
legend(1980,8.0,c("Observed","Fitted"),lty=c(1,2),pch=19,col=c('black','red'))

########
t.c<-seq(from=1,to=142,by=0.01)
T.c<-t.c
C1.c<-cos(2*pi*t.c/12)
S1.c<-sin(2*pi*t.c/12)
C2.c<-cos(2*pi*t.c/6)
S2.c<-sin(2*pi*t.c/6)

y.c<-cbind(1,T.c,C1.c,S1.c,C2.c,S2.c) %*% lWine.coef

plot(t.c,y.c,type="l",ylab=expression(s[t]))
title("Fitted trend and harmonic components")


lWine.res<-lWine-lWine.fitted

lWine.res<-ts(lWine.res, frequency = 12, start = c(1980, 1))
plot(lWine.res,ylab="Monthly Accidental Deaths")
points(lWine.res,pch=19,cex=0.5)
abline(h=0,lty=2)
acf(lWine.res)