library(TSA) data(tempdub) plot(tempdub,type='o') frequency(tempdub) har.=harmonic(tempdub,6) plot(y=har.[,1],x=1:144,type='l') i=0 apply(har.[,1:4],2,function(x){i<<-i+1;lines(y=x,x=1:144,lty=i)}) month.=season(tempdub) # the period sign is included to make the printout from # the commands two line below clearer; ditto below. model2=lm(tempdub~month.-1) # -1 removes the intercept term summary(model2) # Exhibit 3.4 model3=lm(tempdub~month.) # intercept is automatically included so one month (Jan) is dropped summary(model3) Exhibit 3.5 # first creates the first pair of harmonic functions and then fit the model har.=harmonic(tempdub,6) model4=lm(tempdub~har.) summary(model4) har.=harmonic(tempdub,1) model4b=lm(tempdub~har.) summary(model4b) har.=harmonic(tempdub,3) model4a=lm(tempdub~har.) summary(model4a) # Exhibit 3.6 win.graph(width=4.875, height=2.5,pointsize=8) plot(ts(fitted(model4),freq=12,start=c(1964,1)),ylab='Temperature',type='l', ylim=range(c(fitted(model4),tempdub))) # the ylim option ensures that the # y axis has a range that fits the raw data and the fitted values points(tempdub)