ENSO 3.4 Average Data Frame

year <- ENSO[, 1]
ENSO <- c(as.matrix( ENSO[, 2:13]))
y1<- sort(rep(year,12))
# trick to include the leading 0 for months  e.g. 01 
m1 <- rep(1:12,length(year)) + 100
m1 <- substr(format(m1), 2,3)
Edates <- as.Date(paste(y1, m1,"15",sep = "-"))
ENSOtag <- paste(y1, m1, sep = "")
SLOtag <- substr(format(SLO$DATE), 1, 6)

ENSO 3.4 Average Time Series Plot

plot(Edates, ENSO, type="h", xlab = "Date", ylab = "ENSO 3.4 Average", main = "ENSO over Time")

Indicators for Strong/Moderate and El Nino/La Nina

# NAs reported when a match is not possible.
# all places where we have a matching month in SLO with the ENSO subset >=.5
hold1 <- match(SLOtag, ENSOtag[ENSO >= 0.5])
SLOStrongLaNina <- !is.na(hold1)
hold2 <- match(SLOtag, ENSOtag[ENSO <= -0.5])
SLOStrongElNino <- !is.na(hold2)
hold3 <- match(SLOtag, ENSOtag[ENSO > -0.5 & ENSO < 0])
SLOModerateLaNina <- !is.na(hold3)
hold4 <- match(SLOtag, ENSOtag[ENSO < 0.5 & ENSO > 0])
SLOModerateElNino <- !is.na(hold4)
SLOdates <- as.Date(format(SLO$DATE), format = "%Y%m%d")

Precipitation Plotted for Strong and Moderate Storms

set.panel(2,1)
## plot window will lay out plots in a 2 by 1 matrix
#strong El Ninos and La Ninas
#Plotted for the whole time period
plot(SLOdates, SLO[, 2], type="h", xlab = "Date", ylab = "Precipitation (in)", 
     main = "Precipitation over Time")
plot(SLOdates, SLO[, 2], type="n", xlab = "Date", ylab = "Precipitation (in)", 
     main = "Extreme 3.4 Averages")
points(SLOdates[SLOStrongLaNina], SLO[SLOStrongLaNina, 2], col = rgb(0,0,1, alpha = 0.5) , 
       pch = 16)
points(SLOdates, SLO[, 2], type="n")
points(SLOdates[SLOStrongElNino], SLO[SLOStrongElNino, 2], col = rgb(1,0,0, alpha = 0.5) , 
       pch = 16)

#Plotted for a 10 year period to see El Nino and La Nina switch off
trange <- as.Date(c("20000101", "20100101"), format = "%Y%m%d")
plot(SLOdates, SLO[, 2], type="h", xlab = "Date (2000-2010)", ylab = "Precipitation (in)",
     main = "Precipitation over Time", xlim = trange, ylim = c(0.25, 6))
plot(SLOdates, SLO[, 2], type="n", xlab = "Date (2000-2010)", ylab = "Precipitation (in)",
     main = "Extreme 3.4 Averages", xlim = trange, ylim = c(0.25, 6))
points(SLOdates[SLOStrongLaNina], SLO[SLOStrongLaNina, 2], col = rgb(0,0,1, alpha = 0.5) , 
       pch = 16)
points(SLOdates, SLO[, 2], type="n")
points(SLOdates[SLOStrongElNino], SLO[SLOStrongElNino, 2], col = rgb(1,0,0, alpha = 0.5) , 
       pch = 16)

#moderate El Ninos and La Ninas
#Plotted for the whole time period
plot(SLOdates, SLO[, 2], type="h", xlab = "Date", ylab = "Precipitation (in)",
     main = "Precipitation over Time")
plot(SLOdates, SLO[, 2], type="n", xlab = "Date", ylab = "Precipitation (in)",
     main = "Moderate 3.4 Averages")
points(SLOdates, SLO[, 2], type="n")
points(SLOdates[SLOModerateLaNina], SLO[SLOModerateLaNina, 2], col = rgb(0,0,1, alpha = 0.5),
       pch = 16)
points(SLOdates[SLOModerateElNino], SLO[SLOModerateElNino, 2], col = rgb(1,0,0, alpha = 0.5),
       pch = 16)

#Plotted for a 10 year period to see El Nino and La Nina switch off
trange <- as.Date(c("19800101", "19900101"), format = "%Y%m%d")
plot(SLOdates, SLO[, 2], type="h", xlab = "Date (1980-1990)", ylab = "Precipitation (in)",
     main = "Precipitation over Time", xlim = trange, ylim = c(0.25, 6))
plot(SLOdates, SLO[, 2], type="n", xlab = "Date (1980-1990)", ylab = "Precipitation (in)",
     main = "Moderate 3.4 Averages", xlim = trange, ylim = c(0.25, 6))
points(SLOdates, SLO[, 2], type="n")
points(SLOdates[SLOModerateLaNina], SLO[SLOModerateLaNina, 2], col = rgb(0,0,1, alpha = 0.5),
       pch = 16)
points(SLOdates[SLOModerateElNino], SLO[SLOModerateElNino, 2], col = rgb(1,0,0, alpha = 0.5),
       pch = 16)

Histograms with Generalized Pareto Fitted Curves

#Strong NINO
set.panel(2,1)
## plot window will lay out plots in a 2 by 1 matrix
nino <- SLO[SLOStrongElNino, 2]             
nino <- nino[nino > 0]
xgrid <- seq(min(nino), max(nino), length = length(nino))
fit1 <- fevd(nino, threshold = 2, type="GP")
PRCP <- SLO[SLOStrongElNino, 2]
PRCP <- PRCP[PRCP > 0.05]
hist(PRCP, freq = FALSE, main = "Strong El Niños", xlab = "Precipitation (in)", 
     ylab = "Density")
dens <- devd(xgrid, scale = fit1$results$par[1], shape = fit1$results$par[2],
             threshold = fit1$threshold, log = FALSE, type = "GP")
lines(xgrid, dens, col = "red")
#Strong NINA
nina <- SLO[SLOStrongLaNina, 2]
nina <- nina[nina > 0]
xgrid2 <- seq(min(nina), max(nina), length = length(nina))
fit2 <- fevd(nina, threshold = 2, type="GP")
PRCP <- SLO[SLOStrongLaNina, 2]
PRCP <- PRCP[PRCP > 0.05]
hist(PRCP, freq = FALSE, main = "Strong La Niñas", xlab = "Precipitation (in)", 
     ylab = "Density")
dens2 <- devd(xgrid2, scale = fit2$results$par[1], shape = fit2$results$par[2], 
              threshold = fit2$threshold, log = FALSE, type = "GP")
lines(xgrid2, dens2, col = "blue")

#comparing El Nino and La Nina Fitted Curves
set.panel(1,1)
## plot window will lay out plots in a 1 by 1 matrix
plot(-1, -1, xlim = c(0,5), ylim = c(0, 1.8), main = "Strong La Niñas vs. Strong El Niños",
     xlab = "Precipitation (in)", ylab = "Density")
lines(xgrid, dens, col = "red")
lines(xgrid2, dens2, col = "blue")
legend("topright", c("El Niño", "La Niña"), lty = c(1,1), lwd = 1, col = c("red", "blue"))

#Moderate NINO
set.panel(2,1)
## plot window will lay out plots in a 2 by 1 matrix
nino2 <- SLO[SLOModerateElNino, 2]
nino2 <- nino2[nino2 > 0]
xgrid <- seq(min(nino2), max(nino2), length = length(nino2))
fit3 <- fevd(nino2, threshold = 2, type="GP")
PRCP <- SLO[SLOModerateElNino, 2]
PRCP <- PRCP[PRCP > 0.05]
hist(PRCP, freq = FALSE, main = "Moderate El Niños", xlab = "Precipitation (in)", 
     ylab = "Density")
dens <- devd(xgrid, scale = fit3$results$par[1], shape = fit3$results$par[2], 
             threshold = fit3$threshold, log = FALSE, type = "GP")
lines(xgrid, dens, col = "red")
#Moderate NINA
nina2 <- SLO[SLOModerateLaNina, 2]
nina2 <- nina2[nina2 > 0]
xgrid2 <- seq(min(nina2), max(nina2), length = length(nina2))
fit4 <- fevd(nina2, threshold = 2, type="GP")
PRCP <- SLO[SLOModerateLaNina, 2]
PRCP <- PRCP[PRCP > 0.05]
hist(PRCP, freq = FALSE, main = "Moderate La Niñas", xlab = "Precipitation (in)", 
     ylab = "Density")
dens2 <- devd(xgrid2, scale = fit4$results$par[1], shape = fit4$results$par[2], 
              threshold = fit4$threshold, log = FALSE, type = "GP")
lines(xgrid2, dens2, col = "blue")

#comparing El Nino and La Nina Fitted Curves
set.panel(1,1)
## plot window will lay out plots in a 1 by 1 matrix
plot(-1, -1, xlim = c(0,5), ylim = c(0, 1.5), main = "Moderate La Niñas vs. 
     Moderate El Niños", xlab = "Precipitation (in)", ylab = "Density")
lines(xgrid, dens, col = "red")
lines(xgrid2, dens2, col = "blue")
legend("topright", c("El Niño", "La Niña"), lty = c(1,1), lwd = 1, col = c("red", "blue"))

Return Level Plots for El Nino vs. La Nina

#Return Levels
xframe <- matrix(NA, nrow = 500, ncol = 2)          #creating empty matrix
xframe <- as.data.frame(xframe)                     #converting matrix to data frame
yframe <- matrix(NA, nrow = 500, ncol = 2)
yframe <- as.data.frame(yframe)
names(xframe) <- c("Return.Year", "Return.Level")   #assigning names to data frame
names(yframe) <- c("Return.Year", "Return.Level")
x <- SLO[SLOStrongElNino, 2]                        #using strong nino storm indicator
y <- SLO[SLOStrongLaNina, 2]                        #using strong nina storm indicator
pthres <- sum(x > 1)/length(x) 
pthres2 <- sum(y > 1)/length(y) 
for (i in (1:500)){                                #obtaining return levels for 1-500
  xr <- qevd(1 / (pthres*((i)*365)), scale = fit1$results$par[1],  
             shape = fit1$results$par[2], threshold = fit1$threshold,
       type = "GP", lower.tail = FALSE) 
  yr <- qevd(1 / (pthres2*((i)*365)), scale = fit2$results$par[1], 
             shape = fit2$results$par[2], threshold = fit2$threshold,type = "GP", 
             lower.tail = FALSE)
  xframe[i, 1] <- i                               #putting return year in data frame
  xframe[i, 2] <- xr                              #putting return level in data frame
  yframe[i, 1] <- i
  yframe[i, 2] <- yr
}
plot(1, type="n", xlim=c(0, 500), ylim=c(3, 13), main = "Return Levels over Return Years",
     xlab = "Return Year", ylab = "Return Level (precipitation)")
lines(xframe[, 1], xframe[, 2], col = "red", lwd = 1)
lines(yframe[, 1], yframe[, 2], col = "blue", lwd = 1)
legend("bottomright", c("Strong El Niño", "Strong La Niña"), lty = c(1,1), lwd = 1, 
       col = c("red", "blue"))

zframe <- matrix(NA, nrow = 500, ncol = 2)
zframe <- as.data.frame(zframe)
wframe <- matrix(NA, nrow = 500, ncol = 2)
wframe <- as.data.frame(wframe)
names(zframe) <- c("Return.Year", "Return.Level")
names(wframe) <- c("Return.Year", "Return.Level")
z <- SLO[SLOModerateElNino, 2]
w <- SLO[SLOModerateLaNina, 2]
pthres3 <- sum(z > 1)/length(z) 
pthres4 <- sum(w > 1)/length(w) 
for (j in (1:500)){
  zr <- qevd(1 / (pthres3*((j)*365)), scale = fit3$results$par[1],  
             shape = fit3$results$par[2], threshold = fit3$threshold,
             type = "GP", lower.tail = FALSE) 
  wr <- qevd(1 / (pthres4*((j)*365)), scale = fit4$results$par[1], 
             shape = fit4$results$par[2], threshold = fit4$threshold, 
             type = "GP", lower.tail = FALSE)
  zframe[j, 1] <- j
  zframe[j, 2] <- zr
  wframe[j, 1] <- j
  wframe[j, 2] <- wr
}
plot(-1, -1, type="n", xlim=c(0, 500), ylim=c(3, 7),
     main = "Return Levels over Return Years",xlab = "Return Year", 
     ylab = "Return Level (precipitation)")
lines(zframe[, 1], zframe[, 2], col = "red", lwd = 1)
lines(wframe[, 1], wframe[, 2], col = "blue", lwd = 1)
legend("bottomright", c("Moderate El Niño", "Moderate La Niña"), lty = c(1,1), lwd = 1, 
       col = c("red", "blue"))