ANALISI BLOCK MAXIMA PER RX1 STAGIONALI

Serie stagionali ottenute dalle serie di valori giornalieri SCIA.

rm(list = objects())

library("evd")

# 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.delta.lower <- rep(NA, length = length(tempi.ritorno))
    confidence.delta.upper <- rep(NA, length = length(tempi.ritorno))
    confidence.profile.lower <- rep(NA, length = length(tempi.ritorno))
    confidence.profile.upper <- rep(NA, length = length(tempi.ritorno))

    precReturn <- data.frame(idColumn, return.period, return.level, confidence.delta.lower, 
        confidence.delta.upper, confidence.profile.lower, confidence.profile.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)
    index.id <- which(info.anni$idstaz == ID)

    if (!length(index.id)) 
        stop(paste0(ID, "stazione non trovata in file annoI"))
    YYYY <- info.anni[index.id, ]$annoI

    # YYYY: prima di questo anno dobbiamo eliminare la serie
    if (YYYY != 1951) 
        x <- subset(x, x >= YYYY)


    plot(x$yy, x$prcp, type = "l", xlab = "Years", ylab = "mm", main = paste0("Station -- ", 
        ID))
    points(x$yy, x$prcp, pch = 21, bg = "red")

    index <- which(is.na(x$prcp))

    # con pacchetto evd non sarebbe necessario
    if (length(index)) 
        x <- x[-index, ]


    annual.fit <- tryCatch({
        fgev(x$prcp)
    }, 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.delta.lower = NA, confidence.delta.upper = NA, 
            rl.confidence.profile.lower = NA, rl.confidence.profile.upperNA)
        my.shape <- c(NA, NA, NA)
        return(list(precReturn, my.shape))
    })


    plot(annual.fit, which = 1)
    plot(annual.fit, which = 2)
    plot(annual.fit, which = 3)
    plot(annual.fit, which = 4)
    parametri <- round(annual.fit$param, 2)

    cat(sprintf("LOcation Scale and Shape; PACKAGE evd: %s\n", parametri))

    lapply(1:length(tempi.ritorno), function(i) {

        print(sprintf("Calcolo livello di ritorno per tempo: %s anni", tempi.ritorno[i]))
        stima.evd <- fgev(x$prcp, prob = 1/(tempi.ritorno[i]))
        livello.evd <- round(stima.evd$param["quantile"], 2)
        cat(sprintf("RLevel; PACKAGE evd: %s\n", livello.evd))
        profile.object <- evd:::profile.evd(stima.evd, which = "quantile", conf = 0.95)
        plot(profile.object)
        profile.int <- confint(profile.object, level = 0.95, parm = "quantile")
        # intervallo mediante profile likelihood
        profLow <- profile.int[, "lower"]
        profUp <- profile.int[, "upper"]

        # intervalli di confidenza assumendo ipotesi normalità
        profLowW <- confint(stima.evd)[, 1][1]
        profUpW <- confint(stima.evd)[, 2][1]


        precReturn$idColumn[i] <<- ID
        precReturn$return.period[i] <<- tempi.ritorno[i]
        precReturn$return.level[i] <<- livello.evd
        precReturn$confidence.delta.lower[i] <<- profLowW
        precReturn$confidence.delta.upper[i] <<- profUpW

        precReturn$confidence.profile.lower[i] <<- profLow
        precReturn$confidence.profile.upper[i] <<- profUp

    })
    ## fine lapply plot dei tempi di ritorno e intervalli di confidenza (in
    ## rosso) con profile likelihood
    lineaL <- approx(seq(10, 100, 10), precReturn$confidence.profile.lower)
    lineaU <- approx(seq(10, 100, 10), precReturn$confidence.profile.upper)

    plot(annual.fit, which = 4)
    lines(lineaL$x, lineaL$y, col = "red")
    lines(lineaU$x, lineaU$y, col = "red")
    ########################### fine parte su ismev

    # delta intervals for shape
    shapeFi.delta <- confint(annual.fit, parm = "shape")

    # parametro shape e intervalli di confidenza mediante profile likelihood
    profile.object <- profile(annual.fit, which = "shape")
    shapeFi.prof <- confint(profile.object)

    my.shape <- c(ID = as.numeric(ID), xi = parametri[3], xi.delta.dwn = shapeFi.delta[1], 
        xi.delta.up = shapeFi.delta[2], xi.prof.dwn = shapeFi.prof[1], xi.prof.up = shapeFi.prof[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 = "dati1951.csv", sep = ";", head = T, check.names = F)

colonne <- ncol(dati)

risultati.gev <- lapply(2:colonne, analisiGEV)

## Stazione: 6220

plot of chunk unnamed-chunk-1

## LOcation Scale and Shape; PACKAGE evd: 46.26
##  LOcation Scale and Shape; PACKAGE evd: 18.75
##  LOcation Scale and Shape; PACKAGE evd: 0.07
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## RLevel; PACKAGE evd: 91.97
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## RLevel; PACKAGE evd: 108.14
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## RLevel; PACKAGE evd: 117.86
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## RLevel; PACKAGE evd: 124.71
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## RLevel; PACKAGE evd: 129.98
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## RLevel; PACKAGE evd: 133.92
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## RLevel; PACKAGE evd: 133.23
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## RLevel; PACKAGE evd: 142.06
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## RLevel; PACKAGE evd: 144.56
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 100 anni"
## RLevel; PACKAGE evd: 140.43
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "profiling shape"
## Stazione: 6225

plot of chunk unnamed-chunk-1

## LOcation Scale and Shape; PACKAGE evd: 26.84
##  LOcation Scale and Shape; PACKAGE evd: 14.8
##  LOcation Scale and Shape; PACKAGE evd: 0.22
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## RLevel; PACKAGE evd: 69.67
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## RLevel; PACKAGE evd: 88.37
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## RLevel; PACKAGE evd: 99.4
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## RLevel; PACKAGE evd: 108.1
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## RLevel; PACKAGE evd: 116.2
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## RLevel; PACKAGE evd: 123.02
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## RLevel; PACKAGE evd: 126.63
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## RLevel; PACKAGE evd: 128.73
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## RLevel; PACKAGE evd: 135.01
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 100 anni"
## RLevel; PACKAGE evd: 132.67
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "profiling shape"
## Stazione: 6259

plot of chunk unnamed-chunk-1

## LOcation Scale and Shape; PACKAGE evd: 22.12
##  LOcation Scale and Shape; PACKAGE evd: 9.86
##  LOcation Scale and Shape; PACKAGE evd: 0.36
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## RLevel; PACKAGE evd: 56.31
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## RLevel; PACKAGE evd: 74.74
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## RLevel; PACKAGE evd: 86.73
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## RLevel; PACKAGE evd: 96.35
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## RLevel; PACKAGE evd: 103.4
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## RLevel; PACKAGE evd: 108.97
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## RLevel; PACKAGE evd: 114.07
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## RLevel; PACKAGE evd: 125.78
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## RLevel; PACKAGE evd: 120.57
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 100 anni"
## RLevel; PACKAGE evd: 122.47
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "profiling shape"
## Stazione: 6278

plot of chunk unnamed-chunk-1

## LOcation Scale and Shape; PACKAGE evd: 24.81
##  LOcation Scale and Shape; PACKAGE evd: 11.14
##  LOcation Scale and Shape; PACKAGE evd: 0.11
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## RLevel; PACKAGE evd: 53.6
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## RLevel; PACKAGE evd: 63.8
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## RLevel; PACKAGE evd: 70.32
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## RLevel; PACKAGE evd: 75.05
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## RLevel; PACKAGE evd: 78.67
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## RLevel; PACKAGE evd: 82.05
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## RLevel; PACKAGE evd: 84.46
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 80 anni"
## RLevel; PACKAGE evd: 86.8
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 90 anni"
## RLevel; PACKAGE evd: 87.83
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 100 anni"
## RLevel; PACKAGE evd: 91.14
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "profiling shape"
## Stazione: 6429

plot of chunk unnamed-chunk-1

## LOcation Scale and Shape; PACKAGE evd: 33.09
##  LOcation Scale and Shape; PACKAGE evd: 18.43
##  LOcation Scale and Shape; PACKAGE evd: 0.51
## [1] "Calcolo livello di ritorno per tempo: 10 anni"
## RLevel; PACKAGE evd: 111.9
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 20 anni"
## RLevel; PACKAGE evd: 160.23
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 30 anni"
## RLevel; PACKAGE evd: 155.09
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 40 anni"
## RLevel; PACKAGE evd: 178.98
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 50 anni"
## RLevel; PACKAGE evd: 172.55
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 60 anni"
## RLevel; PACKAGE evd: 192.11
## [1] "profiling quantile"

plot of chunk unnamed-chunk-1

## [1] "Calcolo livello di ritorno per tempo: 70 anni"
## RLevel; PACKAGE evd: 256.28
## [1] "profiling quantile"

## Error: need at least two non-NA values to interpolate

plot of chunk unnamed-chunk-1


rl.gevDataFrame <- do.call("rbind", lapply(1:(colonne - 1), function(i) {
    risultati.gev[[i]][[1]]
}))

## Error: object 'risultati.gev' not found

xi.gevDataFrame <- do.call("rbind", lapply(1:(colonne - 1), function(i) {
    risultati.gev[[i]][[2]]
}))

## Error: object 'risultati.gev' not found



write.table(rl.gevDataFrame, file = "rl.gevDataFrame1951.csv", sep = ";", col.names = T, 
    row.names = F, quote = F)

## Error: object 'rl.gevDataFrame' not found

write.table(xi.gevDataFrame, file = "xi.gevDataFrame1951.csv", sep = ";", col.names = T, 
    row.names = F, quote = F)

## Error: object 'xi.gevDataFrame' not found