#R airquality dataset 

#EDA
help(airquality)
summary(airquality)
pairs(airquality, panel = panel.smooth, main = "airquality data")

library(nlme)

fit1 <- gls(Ozone~Wind, data=airquality, na.action=na.omit)

head(airquality)

#Make date variable 
paste(1973, airquality$Month, airquality$Day, sep="-")

airquality$Date <- as.Date(paste(1973, airquality$Month, airquality$Day, sep="-"))

library(lattice)
xyplot(Ozone~Date, airquality)

#or

plot(airquality$Date, airquality$Ozone, main="airquality data")

#Notice temporal trend in data? 

fit2 <- gls(Ozone~Wind*Date, data=airquality, na.action=na.omit)
summary(fit2)

#There is a ACF function in the nlme library, but it does not have 95% CIs, 
#so let's use the time series functions.

par(mfrow=c(2,1))
acf(fit2$residuals)
pacf(fit2$residual)

fit3 <- update(fit2, correlation=corAR1()

AIC(fit2, fit3)

summary(fit3)
par(mfrow=c(1,1)
plot(fit3)

#A few more competitors? What about residuals for fit2? 
plot(fit2)

#Variance increases as a power of the absolute fitted values

fit4 <- update(fit2, weights = varPower())
fit5 <- update(fit3, weights = varPower())