# RPD -- Sec 8.2 -- Quarterly Beer Production

beer <- data.frame(time = 1:32,
                   production = c(36.14, 44.60, 44.15, 35.72,
                                  36.19, 44.63, 46.95, 36.90,
                                  39.66, 49.72, 44.49, 36.54,
                                  41.44, 49.07, 48.98, 39.59,
                                  44.29, 50.09, 48.42, 41.39,
                                  46.11, 53.44, 53.00, 42.52,
                                  44.61, 55.18, 52.24, 41.66,
                                  47.84, 54.27, 52.31, 42.03))


plot(production ~ time, type="b", pch=16, data=beer)  # "b": both points & lines


y <- beer$production

  # Pure annual-cycle second-degree trigonometric model (ignoring trend):

X1 <- with(beer, cbind(1, cos(2*pi*time/4),   sin(2*pi*time/4),
                          cos(2*pi*2*time/4), sin(2*pi*2*time/4)))

print(zapsmall(X1))

X1 <- X1[,-5]  # remove fifth column, which is all zeros

XtX1 <- t(X1)%*%X1

XtXinv1 <- solve(XtX1)

betahat1 <- XtXinv1%*%(t(X1)%*%y)

SSRes1 <- drop(t(y)%*%y - t(betahat1)%*%t(X1)%*%y)

s2.1 <- SSRes1/(dim(X1)[1]-dim(X1)[2])

print(zapsmall(X1))
print(zapsmall(XtX1))
print(betahat1)
print(s2.1)


  # Model augmented with a linear trend:

X2 <- with(beer, cbind(X1,time))

XtX2 <- t(X2)%*%X2

XtXinv2 <- solve(XtX2)

betahat2 <- XtXinv2%*%(t(X2)%*%y)

SSRes2 <- drop(t(y)%*%y - t(betahat2)%*%t(X2)%*%y)

s2.2 <- SSRes2/(dim(X2)[1]-dim(X2)[2])

print(zapsmall(X2))
print(zapsmall(XtX2))
print(betahat2)
print(s2.2)


  # To test whether first-degree model is adequate ...

t0 <- betahat2[4]/sqrt(XtXinv2[4,4]*s2.2)

pval.t0 <- 2*(1-pt(abs(t0),dim(X2)[1]-dim(X2)[2]))

print(t0)
print(pval.t0)


  # To test if any trigonometric terms are needed at all ...

XR <- with(beer, cbind(1,time))

XtXinvR <- solve(t(XR)%*%XR)

betahatR <- XtXinvR%*%(t(XR)%*%y)

SSResR <- drop(t(y)%*%y - t(betahatR)%*%t(XR)%*%y)

F0 <- ((SSResR-SSRes2)/3)/s2.2

pval.F0 <- 1-pf(F0,3,dim(X2)[1]-dim(X2)[2])

print(F0)
print(pval.F0)


  # Compute predicted values for final model:

X <- with(beer, cbind(1, cos(2*pi*time/4), sin(2*pi*time/4), time))

XtXinv <- solve(t(X)%*%X)

betahat <- XtXinv%*%(t(X)%*%y)

yhat <- X%*%betahat


get(getOption("device"))()

matplot(beer$time, cbind(beer$production,yhat), type="b", lty=1, pch=c(16,1),
        col=c("black","green4"))
legend("bottomright", c("Production","Y-hat"), lty=1, pch=c(16,1),
       col=c("black","green4"))