# I 3

day <- as.numeric(format(x$aTime + x$aOffset, "%j", tz="GMT"))
z <- x[day == 11 & x["dest"] == "ORD",]
z <- z[order(z$aTime), ]
lmORD.model <- lm(arrDelay ~ depDelay, data = z)
lmORD.res <- resid(lmORD.model)
# dat junk be nonuniform errorz
plot(z$aTime, lmORD.res^2)
sigmaSq <- mean(lmORD.res^2)
# summary(lmORD.model)$sigma^2
eiej <- sapply(2:nrow(z), function (x) {
  lmORD.res[x] * lmORD.res[x-1]
})
phi <- -1 * log(mean(eiej) / sigmaSq)
# do we need to divide by something to get s^2? or does it divide out...


q <- x[day == 44 & x["dest"] == "ORD",]
q <- q[order(q$aTime), ]
rownames(q) <- NULL

plot(q$depDelay,q$arrDelay)
newLmORD.model <- lm(arrDelay ~ depDelay, data = q)
newLmORD.res <- resid(newLmORD.model)
summary(newLmORD.model)
lm.pred <- predict(newLmORD.model, q, interval="predict")
mean((q$arrDelay > lm.pred[,3]) | (q$arrDelay < lm.pred[,2]))
# [1] 0.05511811

newSigmaSq <- mean(newLmORD.res^2) # divide by something?
covarianceFn <- function(i,j) {newSigmaSq * exp(-1 * phi * abs(i - j))}
Vphi <- outer(1:nrow(q), 1:nrow(q), covarianceFn)
VphiInv <- solve(Vphi)
C <- chol(VphiInv)
q$oldDepDelay <- q$depDelay
q$oldArrDelay <- q$arrDelay
q$depDelay <- C %*% q$depDelay
q$arrDelay <- C %*% q$arrDelay

glmORD.model <- lm(arrDelay ~ depDelay, data = q)
plot(q$depDelay,q$arrDelay)

plot(q$aTime, resid(glmORD.model)^2)

coef(summary(newLmORD.model))[,1]
coef(summary(glmORD.model))[,1]
# differ a lot, but is it because glm is in squiggle space?

plot(q$oldDepDelay, q$oldArrDelay)
glm.pred <- solve(C) %*% predict(glmORD.model, q, interval="predict")
mean((q$arrDelay > glm.pred[,3]) | (q$arrDelay < glm.pred[,2]))