# Studio 12 (Climate Change) ----

# Be sure to read the instructions file!

#-------------------------------
# Problem 1
# See the instructions for this studio.

dataFile <- "./climate.csv"

climateData <- read.csv(dataFile, header = TRUE)
obs <- climateData$obs
ghg <- climateData$ghg
nat <- climateData$nat
aer <- climateData$aer
ant <- ghg + aer
years <- climateData$year

# convert the observations (which are anomalies with respect to the mean of the period 1961-1990) to be anomalies with
# respect to pre-industrial
refp <- 1850:1870
avgRef <- mean(obs[match(refp, years)])
obs <- obs - avgRef

# -------------------------------------
# Problem 1a: Plotting the observations and different forcings
# see the instructions for this studio.
studio12_problem_1a = function() {
  plot(
    years,
    obs,
    col = "black",
    type = "o",
    ylim = c(min(c(obs, ghg, aer, nat)), max(c(obs, ghg, aer, nat))),
    xlab = "Year",
    ylab = "Temperature anomaly [K]"
  )
  lines(years, ghg, col = "red")
  lines(years, aer, col = "blue")
  lines(years, nat, col = "green4")

  legend_labels = c("HadCRUT5", "hist-ghg", "hist-aer", "hist-nat")
  legend(
    "topleft",
    legend_labels,
    col = c("black", "red", "blue", "green4"),
    lty = rep(1, 6),
    pch = c(1, NA, NA, NA, NA, NA)
  )

  # Do not change the above code.
  # ********* YOUR CODE BELOW HERE ***********
  cat("1a: the reason for ... is ...\n")
}

# -------------------------------------
# Problem 1b: Linear regression model 1
# see the instructions for this studio.
studio12_problem_1b = function() {
  # Do not change the above code.
  # ********* YOUR CODE BELOW HERE ***********
  model <- lm(obs ~ ghg + nat + aer)
  print(summary(model))

  cat("problem 1b:\n")
  cat("----------------------------------------\n")
  cat("R^2 = ...\n")
  cat("obs = ...\n")
  cat("least significant is ...\n")
  cat("the reason that ... is the least significant is ...\n")

  # Do not change the line below.
  return(model)
}

# -------------------------------------
# Problem 1c: Linear regression model 2
# see the instructions for this studio.
studio12_problem_1c = function() {
  # Do not change the above code.
  # ********* YOUR CODE BELOW HERE ***********
  model <- lm(obs ~ nat + ant)
  print(summary(model))

  cat("problem 1c:\n")
  cat("----------------------------------------\n")
  cat("R^2 = ...\n")
  cat("obs = ...\n")
  cat("signifcance: ...\n")

  # Do not change the line below.
  return(model)
}

# -------------------------------------
# Problem 1d: Plot fitted values
# see the instructions for this studio.
studio12_problem_1d = function() {
  GhgNatAer <- studio12_problem_1b()
  AntNat <- studio12_problem_1c()

  # Do not change the above code.
  # ********* YOUR CODE BELOW HERE ***********
  fitGhgNatAer <- fitted(GhgNatAer)
  fitAntNat <- fitted(AntNat)

  plot(
    years,
    obs,
    col = "black",
    type = "o",
    ylim = c(min(c(obs, ghg, aer, nat)), max(c(obs, ghg, aer, nat))),
    xlab = "Year",
    ylab = "Temperature anomaly [K]"
  )
  lines(years, ghg, col = "red")
  lines(years, aer, col = "blue")
  lines(years, nat, col = "green4")
  lines(years, fitGhgNatAer, col = "cyan", lwd = 2)
  lines(years, fitAntNat, col = "magenta", lwd = 2)
  leglabs = c("HadCRUT5", "hist-ghg", "hist-aer", "hist-nat", "FittedGhgNatAer", "FittedAntNat")
  legend(
    "topleft",
    leglabs,
    col = c("black", "red", "blue", "green4", "cyan", "magenta"),
    lty = rep(1, 6),
    pch = c(1, NA, NA, NA, NA, NA),
  )
}

# -------------------------------------
# Problem 1e: Attributable warming
# see the instructions for this studio.
studio12_problem_1e = function() {
  model <- studio12_problem_1b()
  bhat <- model$coefficients
  attWarming <- data.frame(ghg = bhat[2] * ghg, nat = bhat[3] * nat, aer = bhat[4] * aer)
  last_row <- tail(attWarming, n=1)

  cat("Attribution of ghg: ", last_row[,'ghg'], "\n")
  cat("Attribution of nat: ", last_row[,'nat'], "\n")
  cat("Attribution of aer: ", last_row[,'aer'], "\n")
}

studio12_problem_1a()
studio12_problem_1b()
studio12_problem_1c()
studio12_problem_1d()
studio12_problem_1e()