Calcolo tempi di ritorno con pacchetto extremes per le serie che iniziano dal 1916 in poi (Serie annuali dei massimi di precipitazione SCIA, rx1)
rm(list = objects())
library("extRemes")
## Loading required package: Lmoments
## Loading required package: distillery
## Loading required package: car
##
## Attaching package: 'extRemes'
##
## The following objects are masked from 'package:stats':
##
## qqnorm, qqplot
# elimina warnings
options(warn = -1)
# funzioni per il calcolo del massimo annuale (RX1) block-maxima approach
analisiGEV <- function(col) {
# data.frame per il livelli di ritorno
idColumn <- rep(NA, length = length(tempi.ritorno))
return.period <- rep(NA, length = length(tempi.ritorno))
return.level <- rep(NA, length = length(tempi.ritorno))
confidence.boot.lower <- rep(NA, length = length(tempi.ritorno))
confidence.boot.upper <- rep(NA, length = length(tempi.ritorno))
precReturn <- data.frame(idColumn, return.period, return.level, confidence.boot.lower,
confidence.boot.upper)
ID <- names(dati)[col]
cat(sprintf("Stazione: %s", ID))
# associamo i dati
x <- dati[, c(1, col)]
names(x) <- c("yy", "prcp")
# subset dei dati sulla base di quanto riportato nel file annoI
file.annoI <- list.files(pattern = "^annoI.+txt$")
if (length(file.annoI) != 1)
stop("Impossibile: ci deve essere un solo file")
# cerchiamo anno inizio
info.anni <- read.csv(file = file.annoI[1], sep = ";", head = T)
names(info.anni)[1] <- "idstaz"
index.id <- which(info.anni$idstaz == ID)
if (!length(index.id)) {
YYYY <- 1951
} else {
YYYY <- info.anni[index.id, ]$annoI
}
# YYYY: prima di questo anno dobbiamo eliminare la serie
if (YYYY == 1951)
return()
x <- subset(x, x >= YYYY)
index <- which(is.na(x$prcp))
# con pacchetto evd non sarebbe necessario
if (length(index))
x <- x[-index, ]
annual.fit <- tryCatch({
fevd(x = prcp, data = x, type = "GEV", location.fun = ~1, use.phi = F,
units = "mm", time.units = "year")
}, error = function(err) {
print(paste(ID, "gev", sep = " "))
print(err)
precReturn <- data.frame(ID = ID, conf.level = NA, return.level = NA,
return.period = NA, confidence.boot.lower = NA, confidence.boot.upper = NA)
my.shape <- c(NA, NA, NA)
return(list(precReturn, my.shape))
})
plot(annual.fit, rperiods = tempi.ritorno)
parametri <- round(annual.fit$results$par, 2)
lapply(1:length(tempi.ritorno), function(i) {
print(sprintf("Calcolo livello di ritorno per tempo: %s anni", tempi.ritorno[i]))
ci.ris <- ci(annual.fit, type = "return.level", return.period = tempi.ritorno[i],
R = 999, alpha = 0.05, method = "boot")
profLow <- round(ci.ris[1], 2)
livello.extr <- round(ci.ris[2], 2)
profUp <- round(ci.ris[3], 2)
precReturn$idColumn[i] <<- ID
precReturn$return.period[i] <<- tempi.ritorno[i]
precReturn$return.level[i] <<- livello.extr
precReturn$confidence.boot.lower[i] <<- profLow
precReturn$confidence.boot.upper[i] <<- profUp
})
ci.shape <- ci(annual.fit, method = "boot", R = 999, type = "parameter",
alpha = 0.05, which.par = 3)
my.shape <- c(ID = as.numeric(ID), xi = round(ci.shape[2], 2), xi.boot.dwn = round(ci.shape[1],
2), xi.boot.up = round(ci.shape[3], 2))
list(precReturn, my.shape)
} #fine analisiGEV
###################### inizio programma
tempi.ritorno <- seq(10, 100, by = 10)
# dati di input: serie che cominciano dal 1951
dati <- read.csv(file = "dati1916.csv", sep = ";", head = T, check.names = F)
colonne <- ncol(dati)
risultati.gev <- lapply(2:colonne, analisiGEV)
## Stazione: 6220
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## [1] "Calcolo livello di ritorno per tempo: 100 anni"
## Stazione: 6278
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## [1] "Calcolo livello di ritorno per tempo: 100 anni"
## Stazione: 6225
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## [1] "Calcolo livello di ritorno per tempo: 100 anni"
rl.gevDataFrame <- do.call("rbind", lapply(1:(colonne - 1), function(i) {
risultati.gev[[i]][[1]]
}))
xi.gevDataFrame <- do.call("rbind", lapply(1:(colonne - 1), function(i) {
risultati.gev[[i]][[2]]
}))
write.table(rl.gevDataFrame, file = "rl.gevDataFrame1951.csv", sep = ";", col.names = T,
row.names = F, quote = F)
write.table(xi.gevDataFrame, file = "xi.gevDataFrame1951.csv", sep = ";", col.names = T,
row.names = F, quote = F)