##### Lesson 60: Extreme Value Distribution in R #####

# Setting the working directory 
setwd(“provide the path to your folder“)

library(locfit)
library(logspline)

## EVD - Simulations ## 

# Normal Origins #
par(mfrow=c(3,1))
x = rnorm(10000,0,1)
plot(logspline(x),xlim=c(-5,5),ylim=c(0,0.5),font=2,cex.lab=2,font.lab=2,xlab="Normal Distribution")

N = 1000
ave = matrix(NA,nrow=N,ncol=1)
max1 = matrix(NA,nrow=N,ncol=1)

for (i in 1:N)
{
  x = rnorm(N,0,1)
  
  lines(locfit(~x),col="grey")
  points(mean(x),0,col="blue",pch=17)
  points(max(x),0,col="red",pch=17)
  ave[i] = mean(x)
  max1[i] = max(x)
  Sys.sleep(0.01)  
}

Sys.sleep(1)
plot(locfit(~ave),xlim=c(-5,5),ylim=c(0,9),ylab="",cex.lab=2,font=2,font.lab=2,xlab="Normal Distribution",col="white")
lines(locfit(~ave),lwd=2,col="blue")

Sys.sleep(1)
plot(locfit(~max1),xlim=c(-5,5),font=2,font.lab=2,ylab="",xlab="Extreme Value Distribution (Gumbel)",cex.lab=2,col="white")
lines(locfit(~max1),lwd=2,col="red")


## EVD from Exponential## 
par(mfrow=c(3,1))
x = rexp(10000,0.5)
hist(x,prob=T,xlim=c(0,25),ylab="",main="",font=2,cex.lab=2,font.lab=2,xlab="Exponential Distribution")
plot(logspline(x),add=T)

N = 1000

ave = matrix(NA,nrow=N,ncol=1)
max1 = matrix(NA,nrow=N,ncol=1)

for (i in 1:N)
{
  x = rexp(N,0.5)
  
  lines(locfit(~x),col="grey")
  points(mean(x),0,col="blue",pch=17)
  points(max(x),0,col="red",pch=17)
  ave[i] = mean(x)
  max1[i] = max(x)
  Sys.sleep(0)  
}

Sys.sleep(1)
plot(locfit(~ave),xlim=c(0,25),ylim=c(0,4),ylab="",cex.lab=2,font=2,font.lab=2,xlab="Normal Distribution",col="white")
lines(locfit(~ave),lwd=2,col="blue")

Sys.sleep(1)
plot(locfit(~max1),xlim=c(0,25),font=2,font.lab=2,ylab="",xlab="Extreme Value Distribution (Gumbel)",cex.lab=2,col="white")
lines(locfit(~max1),lwd=2,col="red")


## EVD from Uniform## 
par(mfrow=c(3,1))
x = runif(100000,0,1)
#hist(x,prob=T,xlim=c(0,1),ylab="",main="",font=2,cex.lab=2,font.lab=2,xlab="Uniform Distribution")
plot(density(x),xlim=c(0,1),ylab="",main="",font=2,cex.lab=2,font.lab=2,xlab="Uniform Distribution")

N = 200

ave = matrix(NA,nrow=N,ncol=1)
max1 = matrix(NA,nrow=N,ncol=1)

for (i in 1:N)
{
  x = runif(N,0,1)
  
  #  lines(locfit(~x),col="grey")
  points(mean(x),0,col="blue",pch=17)
  points(max(x),0,col="red",pch=17)
  ave[i] = mean(x)
  max1[i] = max(x)
  Sys.sleep(0)  
}

Sys.sleep(1)
plot(locfit(~ave),xlim=c(0,1),ylim=c(0,30),ylab="",cex.lab=2,font=2,font.lab=2,xlab="Normal Distribution",col="white")
lines(locfit(~ave),lwd=2,col="blue")

Sys.sleep(1)
plot(locfit(~max1),xlim=c(0,1),ylim=c(0,65),font=2,font.lab=2,ylab="",xlab="Extreme Value Distribution (Weibull)",cex.lab=2,col="white")
lines(density(max1),lwd=2,col="red")


## Plots of Type I, II and III ##

z = sort(runif(1000,-5,5))

alpha = 0
beta = 1
gamma = 1

zprime = (z - alpha)/beta

z1 = zprime
z2 = zprime[zprime>0]
z3 = zprime[zprime<0]

fz1 = exp(-exp(-z1)-z1)

fz2 = gamma*z2^(-gamma-1)*exp(-(1/z2^gamma))

fz3 = (-1/z3)*(gamma)*(-z3)^gamma*exp(-(-z3)^gamma)
  
plot(z1,fz1,type="l",ylim=c(0,1),lwd=2,xlab="z",ylab="f(z)",font=2,font.lab=2)
lines(z2,fz2,col="red",lwd=2)
lines(z3,fz3,col="blue",lwd=2)
abline(v=0)
legend(1, 0.8, legend=c("Gumbel", "Frechet", "Weibull"),
       col=c("black", "red", "blue"), lwd=c(2,2,2))

### GEV and Real Data ###

# fatigue data
fatigue = read.table("cycles_fatigue.txt",header=F)

# temperature data 
temperature_data = read.csv("cp_temperature.csv",header=T)

install.packages("extRemes")
library(extRemes)

#____________#
# Generate Gumbel Random numbers
x = revd(10000,loc=0,scale=1,shape=0)
hist(x,prob=T,xlab="Random Variables from Gumbel (location = 0, scale = 1, shape =0)",main="Gumbel Distribution",ylab="f(x)",font=2,family="serif",font.lab=2,cex.lab=1.5)
plot(logspline(x),add=T)

#____________#
# Generate Frechet Random numbers

# Frechet distribution plot
y = revd(10000,loc=0,scale=1,shape=0.2)
hist(y,prob=T,ylim=c(0,0.4),xlab="Random Variables from Frechet (location = 0, scale = 1, shape = 0.2)",main="Frechet Distribution",ylab="f(x)",font=2,family="serif",font.lab=2,cex.lab=1.5)
plot(logspline(y),add=T,col="red")
plot(logspline(x),add=T)

#____________#
# Generate Frechet Random numbers with large shape parameter

# Frechet distribution plot
y = revd(10000,loc=0,scale=1,shape=0.01)
hist(y,prob=T,ylim=c(0,0.4),xlab="Random Variables from Frechet (location = 0, scale = 1, shape = 0.2)",main="Frechet Distribution",ylab="f(x)",font=2,family="serif",font.lab=2,cex.lab=1.5)
plot(logspline(y),add=T,col="red")
plot(logspline(x),add=T)


#____________#
# Generate Weibull Random numbers with large shape parameter

## Weibull Distribution Plot 
z = revd(10000,loc=0,scale=1,shape=-0.6)
hist(z,prob=T,ylim=c(0,0.5),xlab="Random Variables from Weibull (location = 0, scale = 1, shape = -0.6)",main="Weibull Distribution",ylab="f(x)",font=2,family="serif",font.lab=2,cex.lab=1.5)
plot(logspline(z),add=T,col="red")

# Use FEVD Function #

fatigue_cycles = as.matrix(fatigue)

fit = fevd(fatigue_cycles,type="GEV")
summary(fit)


# Block Maxima Figure #

# Extract 2015 Daily Temperature Data #
index = which((temperature_data[,1] == 2015) | (temperature_data[,1] == 2016) | (temperature_data[,1] == 2017))
temperature_block = jitter(temperature_data[index,4])

plot(temperature_block,xlab="days in 2015, 2016, 2017",ylab="Average Temperature",font=2,font.lab=2, main="Block Maxima (extract the red triangles)")
maxy1 = max(temperature_block[1:365])
maxy2 = max(temperature_block[366:(365+366)])
maxy3 = max(temperature_block[732:1096])

day_max1 = which(temperature_block == maxy1)
day_max2 = which(temperature_block == maxy2)
day_max3 = which(temperature_block == maxy3)

points(day_max1, maxy1,col="red",pch=17,cex=2)
points(day_max2, maxy2,col="red",pch=17,cex=2)
points(day_max3, maxy3,col="red",pch=17,cex=2)

abline(v = c(365,731))
text(200,10,"Year 1",font=2,col="blue")
text(550,10,"Year 2",font=2,col="blue")
text(900,10,"Year 3",font=2,col="blue")

# Annual Maximum Ave Temperature #
yr1 = 1951
yr2 = 2017
n = yr2 - yr1 + 1

annmax = matrix(NA,nrow=n,ncol=1)
for (i in 1:n)
{
  yr = 1950 + i
  index = which((temperature_data[,1] == yr))
  temperature_yr = temperature_data[index,4]
  annmax[i,1] = max(temperature_yr,na.rm=T)
}

hist(annmax,xlab="Annual Maximum Temperature",font=2,main="",font.lab=2)

fit_temperature = fevd(annmax,type="GEV")
plot(fit_temperature)

summary(fit_temperature)