setwd('U:/public_html/588/Class17')    
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)

 anova( lm(temp ~ 1))   # this is "Model 0"

 # this is an "Extra SS" F-test for all 4 terms combined:

 ( F <- (7.93 - 0.98)/4/(0.983/124) )

  prange <- (min(year):2100) - 1940
  new <- data.frame( yr = prange, yr2 = prange^2, sinp = sin(2*pi*prange/pd), cosp = cos(2*pi*prange/pd) )
  
  pred.plim <- predict(lmp, new, interval="prediction")
  pred.clim <- predict(lmp, new, interval="confidence")
  matplot(prange,cbind(pred.clim, pred.plim[,-1]),
        lty=c(1,2,2,3,3), type="l", ylab="predicted y", ylim=c(-0.5,1.5), lwd=2)
    points(yr, temp)