# STA 437/1005 Fall 2010, Assignment 1, Part I.

# Read the data for Sundays into ws and Mondays into wm.

ws <- read.table("ass1-wind-sundays.txt",head=TRUE)
wm <- read.table("ass1-wind-mondays.txt",head=TRUE)


# Compute averages over the 12 locations for Sundays and Mondays.

ws.ave <- rowMeans(ws)
wm.ave <- rowMeans(wm)


# Model the seasonal effect, with no transformation, log, and sqrt.  Also 
# find the fitted values and residuals.  (This is done explictly to show what 
# is happening, but it could be done more easily with the "predict" and
# "residual" functions.)

year.angle <- 2*pi*(((1:nrow(ws))/52.18)%%1)

se.id <- lm(ws.ave ~ sin(year.angle) + cos(year.angle))
print(summary(se.id))
se.fit.id <- coef(se.id)[1] + coef(se.id)[2]*sin(year.angle) +
                              coef(se.id)[3]*cos(year.angle)
se.res.id <- ws.ave - se.fit.id

se.log <- lm(log(ws.ave) ~ sin(year.angle) + cos(year.angle))
print(summary(se.log))
se.fit.log <- coef(se.log)[1] + coef(se.log)[2]*sin(year.angle) +
                                coef(se.log)[3]*cos(year.angle)
se.res.log <- log(ws.ave) - se.fit.log

se.sqrt <- lm(sqrt(ws.ave) ~ sin(year.angle) + cos(year.angle))
print(summary(se.sqrt))
se.fit.sqrt <- coef(se.sqrt)[1] + coef(se.sqrt)[2]*sin(year.angle) +
                                  coef(se.sqrt)[3]*cos(year.angle)
se.res.sqrt <- sqrt(ws.ave) - se.fit.sqrt


# Plot fit and residuals from seasonal model for the three transformations.

pdf("a1-wind-p1.pdf")
par(mfrow=c(2,3))

plot (year.angle, ws.ave, pch=".", main="no transformation")
points (year.angle, se.fit.id, pch=20)
plot (year.angle, log(ws.ave), pch=".", main="log transformation")
points (year.angle, se.fit.log, pch=20)
plot (year.angle, sqrt(ws.ave), pch=".", main="sqrt transformation")
points (year.angle, se.fit.sqrt, pch=20)

plot (year.angle, se.res.id, pch=".", ylab="residuals")
abline(h=0)
plot (year.angle, se.res.log, pch=".", ylab="residuals")
abline(h=0)
plot (year.angle, se.res.sqrt, pch=".", ylab="residuals")
abline(h=0)

dev.off()


# Find variances of residuals in summer - cos(year.angle)<=0 - versus winter.

print(var(se.res.id[cos(year.angle)<=0])/var(se.res.id[cos(year.angle)>0]))
print(var(se.res.log[cos(year.angle)<=0])/var(se.res.log[cos(year.angle)>0]))
print(var(se.res.sqrt[cos(year.angle)<=0])/var(se.res.sqrt[cos(year.angle)>0]))


# Produce histograms and normal QQ plots for seasonal effect residuals with 
# the three transformations.

pdf("a1-wind-p2.pdf")
par(mfrow=c(2,3))

hist(se.res.id, nclass=30, main="no transformation")
hist(se.res.log, nclass=30, main="log transformation")
hist(se.res.sqrt, nclass=30, main="sqrt transformation")

qqnorm(se.res.id)
qqnorm(se.res.log)
qqnorm(se.res.sqrt)

dev.off()

# Set swm.ave to sqrt the average wind speed on Mondays minus seasonal effect.
# This is what we will try to predict.  Set sws.ave to the same for Sundays, 
# and sws to sqrt of individual wind speeds on Sundays minus seasonal effect.

swm.ave <- sqrt(wm.ave) - se.fit.sqrt
sws.ave <- sqrt(ws.ave) - se.fit.sqrt
sws <- sqrt(ws) - se.fit.sqrt


# Look at pairwise scatterplots for sqrt of Sunday wind speeds minus seasonal
# effect, and print the correlation matrix.

pdf("a1-wind-p3.pdf")

plot(sws,pch=".")
print(round(cor(sws),2))

dev.off()


# Find principal components from transformed wind speeds.

pc <- prcomp(sws)
print(summary(pc))
print(round(pc$rotation,2))


# Try predicting average Monday wind speed from average Sunday wind speed.

lm.ave <- lm (swm.ave ~ sws.ave)
print(summary(lm.ave))


# Try predicting average Monday wind speed from all Sunday wind speeds.

lm.all <- lm (swm.ave ~ as.matrix(sws))
print(summary(lm.all))


# Try predicting average Monday wind speed from 1, 2, or 3 PCs.

lm.1pc <- lm (swm.ave ~ pc$x[,1])
print(summary(lm.1pc))

lm.2pc <- lm (swm.ave ~ pc$x[,1:2])
print(summary(lm.2pc))

lm.3pc <- lm (swm.ave ~ pc$x[,1:3])
print(summary(lm.3pc))