#
# This file contains R commands that reproduce the analyses 
# presented in Chapter 9
#
# Read daily temperature data
#
maxtemp<-read.table("../data/maxtemp.data",header=T)
names(maxtemp)<-c("year","month","day","temperature")
n<-dim(maxtemp)[1]
time<-1:n
temperature<-maxtemp$temperature
#
# Time-plot of daily temperature data (Figure 9.1)
#
pdf(file="maxtemp.pdf",height=6,width=8)
plot(time,temperature,type="l",xlab="time (days)",ylim=c(-5,30),cex.axis=1.5,
   cex.lab=1.7)
dev.off()
#
# Weekly black smoke levels in Newcastle upon Tyne, with fitted trend curve (Figure 9.2)
#
smoke<-read.table("../data/blacksmoke.data",header=T)
names(smoke)
pdf(file="blacksmoke.pdf",height=6,width=8)
plot(smoke$week,smoke$BS,type="l",xlab="Week",ylab="Black smoke concentration",log="y",cex.lab=1.5,cex.axis=1.5)
lines(smoke$week,smoke$fit,lty=2,lwd=3)
dev.off()
#
#
# Scatterplots of lagged pairs of daily temperatures (Figure 9.3)
#
plot_range<-range(temperature)
pdf(file="maxtemp_scatter.pdf",paper="special",height=6,width=6)
par(mfrow=c(2,2),pty="s")
for (k in 1:4) {
   x<-temperature[1:(n-k)]
   y<-temperature[(k+1):n]
   plot(x,y,xlab="current",ylab="future",xlim=plot_range,ylim=plot_range,pch=19,cex=0.25)
   lines(plot_range,plot_range)
   }
dev.off()
#
# Fit simple harmonic regressionm model to describe seasonal trend
#
xc<-cos(2*pi*time/366)
xs<-sin(2*pi*time/366)
fit<-lm(temperature~xc+xs)
summary(fit)
#
# Calculate residuals from fitted seasonal trend and make scatterplots 
# of lagged pairs of residuals (Figure 9.4)
#
residuals<-fit$resid
n<-length(residuals)
plot_range<-range(residuals)
pdf(file="maxtemp_residuals_scatter.pdf",height=6,width=8)
par(mfrow=c(2,3),pty="s")
for (k in 1:3) {
   x<-residuals[1:(n-k)]
   y<-residuals[(k+1):n]
   plot(x,y,xlab="current",ylab="future",xlim=plot_range,ylim=plot_range,pch=19,cex=0.25,cex.axis=1.5)
   lines(plot_range,plot_range)
   }
for (k in c(7,14,28)) {
   x<-residuals[1:(n-k)]
   y<-residuals[(k+1):n]
   plot(x,y,xlab="current",ylab="future",xlim=plot_range,ylim=plot_range,pch=19,cex=0.25,cex.axis=1.5)
   lines(plot_range,plot_range)
   }
dev.off()
#
# Correlograms of daily temperatures and of residuals from simple harmonic regression model 
# (Figure 9.5)
#
pdf(file="maxtemp_correlograms.pdf",height=6,width=8)
par(mfrow=c(1,2))
r1<-acf(temperature,lag.max=28,plot=F)
plot(r1$lag,r1$acf,pch=19,cex.axis=1.5,xlab="lag",ylab="autocorrelation",ylim=c(-0.2,1))
lines(r1$lag,r1$acf)
lines(c(0,28),c(1,1)*2/sqrt(366),lty=2)
lines(c(0,28),-c(1,1)*2/sqrt(366),lty=2)
r2<-acf(residuals,lag.max=28,plot=F)
plot(r2$lag,r2$acf,pch=19,cex.axis=1.5,xlab="lag",ylab="autocorrelation",ylim=c(-0.2,1))
lines(r2$lag,r2$acf)
lines(c(0,28),c(1,1)*2/sqrt(366),lty=2)
lines(c(0,28),-c(1,1)*2/sqrt(366),lty=2)
dev.off()
#
# Daily temperature data with fitted seasonal trend (Figure 9.6)
#
fitted<-temperature-residuals
pdf(file="maxtemp_fit.pdf",height=6,width=8)
plot(time,temperature,type="l",xlab="time (days)",ylim=c(-5,30),cex.axis=1.5,
   cex.lab=1.7)
lines(time,fitted,lty=2,lwd=2)
dev.off()
#
# Comparison between observed and fitted autocorrelation functions (Figure 9.7)
#
pdf(file="fitted_correlogram.pdf",height=6,width=8)
r2<-acf(residuals,lag.max=28,plot=F)
plot(r2$lag,r2$acf,pch=19,cex.axis=1.5,cex.lab=1.5,xlab="lag",ylab="autocorrelation",ylim=c(-0.2,1))
lines(r2$lag,r2$acf)
lines(0:28,0.684^(0:28),lty=2,lwd=2)
lines(c(0,28),c(1,1)*2/sqrt(366),lty=2)
lines(c(0,28),-c(1,1)*2/sqrt(366),lty=2)
dev.off()
#
# Mean square precdiction errors when forecasting daily temperatures (Section 9.5)
#
predictions<-function(k,rho=0.683) {
   observed<-temperature[(k+1):n]
   predicted<-fitted[(k+1):n]+(rho^k)*residuals[1:(n-k)]
   MSE<-sum((observed-predicted)^2)/(n-k)
   list(predicted=predicted,MSE=MSE)
   }
MSE.results<-matrix(0,3,7)
kvec<-c(1,2,3,7,14,28)
MSE.results[,1]<-1:3
for (j in 1:6) {
   k<-kvec[j]
   MSE.results[1,j+1]<-predictions(k,0.683)$MSE
   MSE.results[2,j+1]<-predictions(k,0)$MSE
   MSE.results[3,j+1]<-predictions(k,1)$MSE
   }
round(MSE.results,2)
#
#     [,1] [,2] [,3]  [,4]  [,5]  [,6]  [,7]
#[1,]    1 4.03 5.96  6.73  7.49  7.75  7.86
#[2,]    2 7.56 7.57  7.58  7.60  7.73  7.86
#[3,]    3 4.78 8.15 10.07 13.20 17.86 16.44
# 
# A simulated time series with a complex seasonal trend, and apprpocximations
# obtained by fitting different numbers of sine-cosine waves (Figures 9.8, 9.9, 9.10)
# 
n<-250
sigma<-0.1
cosmat<-matrix(0,10,250)
ampvec<-c(0.12,0.01,0.2,0.1,0.25,0.05,0.12,0.06,0.03,0.04)
phasevec<-c(0,1,0.5,1.5,0.2,2.0,1.3,0.7,0.9,0.25)
for (i in 1:10) {
   cosmat[i,]<-ampvec[i]*cos((2*pi*i*(1:n)/n)+phasevec[i])
   }
y<-apply(cosmat,2,sum)+sigma*rnorm(n)
c1<-cos(2*pi*(1:n)/n)
s1<-sin(2*pi*(1:n)/n)
c2<-cos(2*pi*2*(1:n)/n)
s2<-sin(2*pi*2*(1:n)/n)
f2<-lm(y~c1+s1+c2+s2)$fitted.values
c3<-cos(2*pi*3*(1:n)/n)
s3<-sin(2*pi*3*(1:n)/n)
c4<-cos(2*pi*4*(1:n)/n)
s4<-sin(2*pi*4*(1:n)/n)
f4<-lm(y~c1+s1+c2+s2+c3+s3+c4+s4)$fitted.values
c5<-cos(2*pi*5*(1:n)/n)
s5<-sin(2*pi*5*(1:n)/n)
c6<-cos(2*pi*6*(1:n)/n)
s6<-sin(2*pi*6*(1:n)/n)
c7<-cos(2*pi*7*(1:n)/n)
s7<-sin(2*pi*7*(1:n)/n)
c8<-cos(2*pi*8*(1:n)/n)
s8<-sin(2*pi*8*(1:n)/n)
f8<-lm(y~c1+s1+c2+s2+c3+s3+c4+s4+c5+s5+c6+s6+c7+s7+c8+s8)$fitted.values
pdf(file="spectral_approx_1.pdf",paper="special",height=4,width=6)
plot(1:n,y,type="l",xlab="time",ylab="series")
lines(1:n,f2,lty=2,lwd=2)
dev.off()
pdf(file="spectral_approx_2.pdf",paper="special",height=4,width=6)
plot(1:n,y,type="l",xlab="time",ylab="series")
lines(1:n,f4,lty=2,lwd=2)
dev.off()
pdf(file="spectral_approx_3.pdf",paper="special",height=4,width=6)
plot(1:n,y,type="l",xlab="time",ylab="series")
lines(1:n,f8,lty=2,lwd=2)
dev.off()