## this is the re-analysis of the Global temperature data, but taking into account the correlated rediuals
#

setwd('C:/local/latest/classes/588/')
data <- read.table("global-anom.txt")
plot(data)

 year <- data[,1]
 temp <- data[,2]

 yr <- year - 1940    # "coded" data
 yr2 <- yr^2          
 
 # analyzing the 64-year period model
  pd <- 64
   sinp <- sin(2*pi*yr/pd)                           
    cosp <- cos(2*pi*yr/pd)     
   lmp <- lm(temp ~ yr + yr2 + sinp + cosp)     
  summary(lmp)

  plot(yr, temp)
    lines(yr, lmp$fit, col="blue", lwd=2)

 # learning the time series structure 

 resid <- lmp$res

 par(mfrow=c(1,2)) 
  acf(resid)
  pacf(resid)   # useful to determine the order of AR model to fit
  
 n <- length(resid)
 r <- cor(resid[1:(n-1)], resid[2:n])              # shortcut  cor(resid[-1],resid[-n])
 
 # will the fit change if we use GLS? 

 rowr <- 0:(n-1)             # powers of r in the first row
 Sigmat <- toeplitz( exp(rowr*log(r)))   # Sigma-matrix in correlation form
 round(Sigmat[1:5,1:5],3)   
 Sinv <- solve(Sigmat)

 X <- cbind(rep(1,n), yr, yr2, sinp, cosp)   # X-matrix
 XtSXi <- solve(t(X) %*% Sinv %*% X)                
 (betaG <- XtSXi %*% t(X) %*% Sinv %*% temp)
 predG <- X %*% betaG

  plot(yr, temp)
    lines(yr, lmp$fit, col="blue", lwd=2)
    lines(yr, predG, col="red", lwd=2)        # not much of a difference!