Esercizio_4

#   ESERCIZIO 4

rm(list = ls())
library(moments)
library(evd)
library(goftest)

data <- read.table("Esercizio4.txt", sep = '\t', header = T)
df <- data[ ,2:6]
year <- data$Anno

#Vettore di durate
duration <- c(1, 3, 6, 12, 24)
par <- matrix(NA, nrow = length(duration), ncol = 2)
colnames(par) <- c("scale", "location")

#METODO DEI QUANTILI REGOLARIZZATI ----------------

for (i in 1:length(duration)) {
  x <- df[ ,i]
  mu <- mean(x)
  sd <- sqrt(var(x))
  scale <- sqrt(6) * sd / pi
  loc <- mu - 0.5772 * scale
  par[i, 1] <- scale
  par[i, 2] <- loc
}

PR <- c(10, 25, 50, 100)
FR <- (1 - 1/PR)
LSPP <- matrix(NA, nrow = length(duration), ncol = length(PR))
for (i in 1:length(duration)) {
  qg <- qgumbel(FR, loc = par[i, 2], scale = par[i, 1])
  LSPP[i, ] <- qg
}

plot(duration, LSPP[ ,1], ylim = c(40, 310), xlab = "Duration [h]", ylab = 'Quantile', main = 'Quantili regolarizzati')
points(duration, LSPP[ ,2])
points(duration, LSPP[ ,3])
points(duration, LSPP[ ,4])

# per ogni tempo di ritorno fitto le rette nel piano logaritmico e poi ritrasformo i parametri ottenuti
par2 <- matrix(NA, nrow = length(PR), ncol = 2)
colnames(par2) <- c("a", "n")
for (i in 1:length(PR)) {
  y <- log(LSPP[ ,i])
  x <- log(duration)
  fit <- lm(y ~ x)
  a <- exp(fit$coefficients[1])
  n <- fit$coefficients[2]
  par2[i, 1] <- a
  par2[i, 2] <- n
}

dd <- 1:24
plot(dd, par2[1, 1] * dd ^ par2[1, 2], ylim = c(40, 310), ylab="depth (mm)", xlab="duration (h)", main="DDF-Regularized quantile method", col = 'red', lwd = 3, type = 'l')
lines(dd, par2[2, 1] * dd ^ par2[2, 2],col = 'blue', lwd = 3)
lines(dd, par2[3, 1] * dd ^ par2[3, 2],col = 'orange', lwd = 3)
lines(dd, par2[4, 1] * dd ^ par2[4, 2],col = 'green', lwd = 3)
legend("topleft", legend=c("Return period=10", "Return period=25", "Return period=50", "Return period=100"),
       col=c("red", "blue", "orange", "green"), lwd=c(1.5, 1.5, 1.5, 1.5), text.width = 0.4, bty = 'n', text.font=3)



# METODO SCALA INVARIANTE -------------

mean_duration <- apply(df, 2, mean)
X <- c()

for(i in 1: length(duration)){
  X <- c(X, df[,i]/mean_duration[i])
}

# APPLICAZIONE L-MOMENTS PER STIMARE PARAMETRI GEV

xs <- sort(X)
N <- length(X)
rank <- seq(1, N)
PP <- rank / (N + 1)


M0 <- 1
xPP <- xs * PP
M1 <- (1/N) * sum(xPP)
xPP2 <- xs * PP^2
M2 <- (1/N) * sum(xPP2)

L1 <- M0
L2 <- 2 * M1 - M0
L3 <- 6 * M2 - 6 * M1 + M0

c <- 2 * L2 / (L3 + 3 * L2) - log(2) / log(3)
shape <- 7.8590 * c + 2.955 * c ^ 2
scale2 <- shape * L2 / ((1 - 2 ^ (-shape)) * gamma(1 + shape))
loc2 <- L1 - scale2 / shape * (1 - gamma(1 + shape))

x2 <- log(duration)
y2 <- log(mean_duration)
fit <- lm(y2 ~ x2)
a2 <- exp(fit$coefficients[1])
n2 <- fit$coefficients[2]

LSPP2 <- matrix(NA, nrow=length(PR), ncol=length(duration))
for(i in 1: length(PR)){
  xT <- qgev(FR[i], scale=scale2, loc=loc2, shape=shape)
  LSPP2[i,] <- xT
}


plot(dd, LSPP2[1, 1] * a2 * dd ^ n2, type="l", ylab="depth (mm)", xlab="duration (h)", main="DDF-scale invariant method",
     col="red", ylim=c(40,310), lwd = 3)
lines(dd, LSPP2[2, 1]* a2 * dd ^ n2, type="l", col="blue", lwd = 3)
lines(dd, LSPP2[3, 1]* a2 * dd ^ n2, type="l", col="orange", lwd = 3)
lines(dd, LSPP2[4, 1]* a2 * dd ^ n2, type="l", col="green", lwd = 3)
legend("topleft", legend=c("Return period=10", "Return period=25", "Return period=50", "Return period=100"),
       col=c("red", "blue", "orange", "green"), lwd=c(1.5, 1.5, 1.5, 1.5), text.font=3, text.width = 0.4, bty = 'n')