mw <- read.csv("http://stat.ufl.edu/~winner/data/miami_temp_rain_month.csv")
attach(mw); names(mw)

## Y = [meantemp , totprcp], center series month
Y <- cbind(meantemp,totprcp)
ser.mnth <- ser.mnth - mean(ser.mnth)

## Fit Regression of Y on Seasonal predictors, ser.mnth, ser.mnth^2 trend
mw.mod1 <- lm(Y ~ Winter + Spring + Summer + ser.mnth + I(ser.mnth^2))
summary(mw.mod1)
summary(manova(mw.mod1), test="Wilks")

E <- resid(mw.mod1)
(SSCP.ERR <- t(E) %*% E)

n <- nrow(E)
(Sigma.hat <- (1/n) * SSCP.ERR)

## Test whether the Linear and quadratic series month trend terms are all 0

mw.mod2 <- lm(Y ~ Winter + Spring + Summer)
summary(mw.mod2)
summary(manova(mw.mod2), test="Wilks")

E1 <- resid(mw.mod2)
(SSCP.ERR1 <- t(E1) %*% E1)

n <- nrow(E1)
(Sigma.hat1 <- (1/n) * SSCP.ERR1)

anova(mw.mod2, mw.mod1, test="Wilks")

par(mfrow=c(2,1))
plot(E1[,1] ~ ser.mnth, pch=16)
abline(lm(E1[,1] ~ ser.mnth),col="red",lwd=4)
plot(E1[,2] ~ ser.mnth, pch=16)
abline(lm(E1[,2] ~ ser.mnth),col="red",lwd=4)


###################################################################
### Matrix Form                                                ####
###################################################################

Z1 <- cbind(rep(1,n),Winter,Spring,Summer)
Z2 <- cbind(ser.mnth, (ser.mnth^2))
Z <- cbind(Z1,Z2)

## Full Model 

(beta.hat <- solve(t(Z) %*% Z) %*% t(Z) %*% Y)
Y.hat <- Z %*% beta.hat
E.m <- Y - Y.hat
(SSCP.ERR.m <- t(E.m) %*% E.m)
(Sigma.hat.m <- (1/n) * SSCP.ERR.m)

## Reduced Model

(beta1.hat <- solve(t(Z1) %*% Z1) %*% t(Z1) %*% Y)
Y1.hat <- Z1 %*% beta1.hat
E1.m <- Y - Y1.hat
(SSCP.ERR1.m <- t(E1.m) %*% E1.m)
(Sigma.hat1.m <- (1/n) * SSCP.ERR1.m)

## F-test based on Wilks' Lambda   (H0: Beta2=0)
### L_{q x (r+1)} = [0,0,0,0,1,0  // 0,0,0,0,0,1]
r <- ncol(Z) - 1
q <- ncol(Z1) - 1 
m <- ncol(Y)

(H.test <- SSCP.ERR1.m - SSCP.ERR.m)
(E.test <- SSCP.ERR.m)
(HE.test <- H.test + E.test)

L <- matrix(c(0,0,0,0,1,0, 0,0,0,0,0,1), byrow=T, ncol=6)

(q1 <- qr(HE.test)$rank)
(q2 <- qr((L %*% solve(t(Z) %*% Z) %*% t(L)))$rank)
(nu <- n - (r + 1))
(s <- min(q1,q2))
(m1 <- 0.5 * (abs(q1 - q2) - 1))
(m2 <- 0.5 * (nu - q -1))
(lambda <- nu - 0.5 * (q1 - q2 + 1))
(u <- 0.25 * (q1*q2 - 2))
t1 <- q1^2*q2^2 - 4
t2 <- q1^2 + q2^2 - 5
(t <- ifelse(t2>0, sqrt(t1/t2), 1))

Wilks.L <- det(E.test) / det(HE.test)
F.W <- ((1 - Wilks.L^(1/t)) / (Wilks.L^(1/t))) * ((lambda*t - 2*u) / (q1*q2))
df1.W <- q1*q2
df2.W <- lambda*t - 2*u 
CV.W <- qf(.95, df1.W, df2.W)
PV.W <- 1 - pf(F.W, df1.W, df2.W)

Wilks.out <- cbind(Wilks.L, df1.W, df2.W, F.W, CV.W, PV.W)
colnames(Wilks.out) <- c("Lambda", "df", "df", "F", "F(.05)", "P-value")
round(Wilks.out, 3)