library(datasets)

plot(USAccDeaths,ylab="Monthly Accidental Deaths")
points(USAccDeaths,pch=19)
title("Accidental Deaths in the US: Monthly totals \n Jan 1973 - Dec 1978")

mean(USAccDeaths)


Xt<-diff(diff(USAccDeaths),lag=12)
plot(Xt,ylab="Differenced Monthly Accidental Deaths")
points(Xt,pch=19)

t.v<-c(1:72)
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)

plot(t.v,C1,type="l")
lines(t.v,S1)
lines(t.v,C2)
lines(t.v,S2)

USA.ACC<-data.frame(as.numeric(USAccDeaths),S1,C1,S2,C2)

par(mfrow=c(2,1))
USA.ACC.fit<-lm(USAccDeaths~S1+C1+S2+C2,data=USA.ACC)$fitted
USA.ACC.coeff<-lm(USAccDeaths~S1+C1+S2+C2,data=USA.ACC)$coeff

plot(USAccDeaths,ylab="Monthly Accidental Deaths",ylim=range(6500,12000))
points(USAccDeaths)
title("Accidental Deaths in the US: Monthly totals \n Jan 1973 - Dec 1978")
USAFitDeaths<-ts(USA.ACC.fit, frequency = 12, start = c(1973, 1))
points(USAFitDeaths,pch=3)
lines(USAFitDeaths,lty=2)
legend(1974,12000,c("Observed","Fitted"),lty=c(1,2))


t.c<-seq(from=1,to=72,by=0.01)
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<-USA.ACC.coeff[1]+USA.ACC.coeff[2]*S1.c+USA.ACC.coeff[3]*C1.c+USA.ACC.coeff[4]*S2.c+USA.ACC.coeff[5]*C2.c

plot(t.c,y.c,type="l",ylim=range(6500,12000),ylab=expression(s[t]))
title("Fitted harmonic components")

par(mfrow=c(1,1))

USA.resid<-USAccDeaths-USAFitDeaths
acf(USA.resid)

USFitDeaths.res<-ts(USA.resid, frequency = 12, start = c(1973, 1))
plot(USFitDeaths.res,ylab="Monthly Accidental Deaths")
points(USFitDeaths.res,pch=19)
abline(h=0,lty=2)