setwd('U:/public_html/588/Class25/')    
X <- read.csv('TucsonAZ.csv', header = T)
X <- X[(X$year >= 1998) & (X$year <= 2003) ,]   

attach(X)
plot(temp, type="l")
n <- dim(X)[1]
t <- 1:n
sint <- sin(2*pi*t/365.25)
cost <- cos(2*pi*t/365.25)  

lm1 <- lm(temp ~ t + sint + cost)
summary(lm1)

## illustrating autocorrelation
 resid <- lm1$resid
 plot(resid[1:(n-1)], resid[2:n] , xlab="Current", ylab= "Next") 

r <-  cor(resid[1:(n-1)], resid[2:n])   # lag 1 autocorrelation

# now doing GLS with matrix Sigmat

nn <- 5  

 rowr <- 0:(nn-1)             # powers of r in the first row
 Sigmat <- toeplitz( r^rowr) # Sigma-matrix in correlation form for AR(1) model
       # we'll find var(Xbar) under assumption of AR(1) model
Sinv <-  solve(Sigmat)
round(Sinv[1:5,1:5],3)

## let's "extrapolate" the answer 
 rowr <- 0:(n-1)             # powers of r in the first row
 Sigmat <- toeplitz( r^rowr)  
 In1 <- diag(n+1)
 onebelow <- In1[2:(n+1),1:n]
 Sinv <- diag(n)*3.461 + onebelow*(-1.657) + t(onebelow)*(-1.657) 
 Sinv[1,1] <- 2.230
 Sinv[n,n] <- 2.230

Xmat <- cbind(1, t, sint, cost) 
invXSX <- solve(t(Xmat) %*% Sinv %*% Xmat)

betahat <- invXSX %*% t(Xmat) %*% Sinv %*% temp
round(betahat,4)
 s <- sqrt(sum( (temp - Xmat %*% betahat)^2)/(n-4))
 Var <- s^2* invXSX
 SE <- sqrt(diag(Var))   ## these are SE(beta); compare this to the lm1 output
 round(SE,5)
  Tstat <- betahat/SE 
  round(Tstat,2)
                        
(Ymean <- tapply(temp, year, mean) )