# Time Series Regression of RSS Temp Anomaly Trend, SATO Index, NINO3.4 Index
## Monthly anomaly data
## D Kelly O'Day; http://processtrends.com; http://chartsgraphs.wordpress.com
## Feb 11-13, 2009
################################################################################
## 1: Setup
################################################################################
link <- "http://processtrends.com/R_Files/Consol_GISS_SATO_Nino34A_RSS.csv"
   options(digits=6)
  par(ps=9); par(cex.axis=c(0.85)); par(cex.lab=c(0.9));par(las="1")
  par(mfrow=c(2,1))
  par(oma=c(3,3,4,2));  par(mar=c(0,5,0,0)); par(bty="n")
  t_cex = 0.9                 # control plot text note font size
################################################################################
## 2: Read Data
################################################################################
raw_data <- read.csv(link,
             skip = 1, sep = ",", dec=".",
             row.names = NULL, header = FALSE,
             colClasses = rep("numeric",7),
             comment.char = "#", na.strings = c("*", "-",-99.9, -999.9),
             col.names = c("Yr", "Mo","Yr_frac", "GISS", "SATO", "NINO34", "RSS") )
in_data <- subset(raw_data, Yr >= 1979)
################################################################################
## 3: Create time series variables
################################################################################
  # Enter lag months for Nino34 and SATA
    n_lag <- -5
    s_lag <- -3
  # Try several lags to find maximum r2
 ts_rss <- ts(RSS, start =c(1979,1), frequency = c(12))
 ts_nino34 <- ts(NINO34,start = c(1979,1), frequency = c(12))
 ts_sato <- ts(SATO, start = c(1979, 1), frequency = c(12))
 ts_yr_frac <- ts(Yr_frac, start = c(1979,1), frequency = 12)
 ts_data <- data.frame(ts_rss, ts_nino34, ts_sato,ts_yr_frac)
  nino_l <- lag(ts_nino34, n_lag)
  sato_l <- lag(ts_sato, s_lag)
m_d <-ts.union(ts_yr_frac, ts_rss, ts_nino34,nino_l, sato_l)
################################################################################
## 4: Develop simple trend line regression model
################################################################################
   trend_lm_fit <- lm(ts_rss ~ ts_yr_frac, data = m_d)
   trend_sum <- summary(trend_lm_fit)
   # get tend line coefficients
   trend_a <- coef(trend_lm_fit)[1]
   trend_b <- coef(trend_lm_fit)[2]
   trend_r2 <-  trend_sum$adj.r.squared
################################################################################
## 5: Multiple Regression: RSS = a + b*yr_frac + c*NINO34_lag + d*SATO_lag
################################################################################
 lm_fit <- lm(ts_rss ~ ts_yr_frac + nino_l + sato_l, data = m_d)
 reg_sum <- summary(lm_fit)
 # get regression coefficients
    a <- coef(lm_fit)[1]
    b <- coef(lm_fit)[2]
    c <- coef(lm_fit)[3]
    d <- coef(lm_fit)[4]
    r2 <-  reg_sum$adj.r.squared
# produce multiple regression estimate
  mr_est <- a +b*ts_yr_frac + c*nino_l + d* sato_l
################################################################################
## 6: Plot 1: Simple Trend Line
################################################################################
 trend_f = paste("RSS = ", signif(trend_a, 3), " +", signif(trend_b,3), "*Yr_Frac")
 plot(ts_rss, col = "grey", xlim = c(1979, 2010),
        ylim = c(-1,1), xaxs="i", yaxs ="i",
        ylab = expression(paste("RSS Anomaly - ", degree, "C")),
        axes = FALSE)
   axis(1, at= FALSE,  labels=FALSE, col = "grey")
    axis(2, labels = TRUE, col = "grey")
     axis(3, labels=FALSE, at = FALSE, col = "grey")
  abline(trend_lm_fit, col = "red")
   abline(v = c(1985, 1990, 1995, 2000,2005, 2010), col = "lightgrey")
    abline(h=1, col = "grey")
    abline(h=-1, col = "grey")
  rect(2001.,-0.5, 2006, -0.25, border = NA, col = "white")
    text(2008,-0.3, paste("adj r squared = ", signif(trend_r2,3)), cex = t_cex, pos = 2)
  rect( 1980, 0.7, 1989.5, 0.9, bor = NA, col = "white")
   text(1979, 0.8, "a) Simple Trend Line", col = "black", pos = 4,adj = 0)
  rect(1979.2, -0.9, 1992, -0.7, border = NA, col = "white")
   text(1982, -0.8, trend_f, col = "black", pos = 4, cex = 0.85, adj = 0)
  lines(c(1980, 1981.7), c(-0.8, -0.8), col = "red")
################################################################################
## 7: PLOT 2:  Multiple Regression
################################################################################
 par(mar=c(1,5,0,0))
 rss_f = paste("RSS = ", signif(a, 3), " +", signif(b,3), "*Yr_Frac + ",signif(c,3), "*NINO34 ", signif(d,3), "*SATO")
  plot(ts_rss, col = "grey", xlim = c(1979, 2010),
     ylim = c(-1, 1),  xaxs="i", yaxs ="i", xlab = "",
     ylab = expression(paste("RSS Anomaly - ", degree, "C")),
     mgp = c(3,1,0))
   abline(h=-1, col = "grey")
    axis(1, col = "grey", at= NULL, labels = TRUE, cex = 1)
     axis(2, col = "grey")
      abline(v = c(1985, 1990, 1995, 2000,2005, 2010), col = "lightgrey")
  points(mr_est, type ="l", col = "red")
  rect( 1980, 0.7, 1997, 0.9, bor = NA, col = "white")
   text(1979, 0.8, "b) Multiple Regression (Yr, NINO34, SATO)", pos = 4, adj = 0)
  rect(2000.2, -0.6, 2008, -0.2, border = NA, col = "white")
   text(2008,-0.3, paste("adj r squared = ", signif(r2,3)), cex = t_cex, pos = 2)
  text(2008,-0.4, paste("Nino34 lag ", n_lag, " months"), cex = t_cex, pos = 2)
   text(2008,-0.5, paste("SATO lag ", s_lag, " months"), cex = t_cex, pos = 2)
  rect(1979.2, -0.9, 2004, -0.7, border = NA, col = "white")
   text(1982, -0.8, rss_f, col = "black", pos = 4, adj = 0, cex = 0.85)
  lines(c(1980, 1981.7), c(-0.8, -0.8), col = "red")
#############################################################################
## 8: Overall Chart Title, footer ##############################################
################################################################################
# Calculate mo, yr & value for last readings
      nr <- nrow(ts_data)                    # Find number of data rows
      mo <- Mo[nr]                           # Find last month of data
      yr <- Yr[nr]                           # Find last year of data
      mo_names <- c("Jan", "Feb", "Mar", "Apr", "May", "June", "July", "Aug", "Sept",  "Oct", "Nov", "Dec")
# Generate and add Title_1, Title_2
    Title_1 <- "RSS Temperature Anomalies: a) Simple Trend Line, b) Multipe Regression"
    Title_2 <- paste("Monthly Data: Jan, 1979 through ", mo_names[mo], ",", yr)
   mtext(Title_1, 3,2, outer = TRUE)
   mtext(Title_2, 3,1, outer = TRUE, cex = 0.85)
# Generate and add bottom footer KOD & System date notes
       mtext("D Kelly O'Day - http://chartgraphs.wordpress.com", 1,1, adj = 0, cex = 0.8, outer=TRUE)
       mtext(format(Sys.time(), "%m/%d/ %Y"), 1, 1, adj = 1, cex = 0.8, outer=TRUE)
#################################################################################