Verifica presenza trend

Analisi del trend nelle serie annuali dei massimi di precipitazione (RX1).

I trend vengono valutati mediante lo stimatore di TheilSenn e assumendo un modello GEV non stazionario.

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

library("zyp")

## Loading required package: Kendall


## elimina warnings
options(warn = -1)


analisiGEV <- function(col) {

    par(oldPar)

    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)) {
        YYYY <- 1951
    } else {
        YYYY <- info.anni[index.id, ]$annoI
    }

    ## YYYY: prima di questo anno dobbiamo eliminare la serie Consideriamo solo
    ## le serie che partono dal 1951
    if (YYYY != 1951) {
        sprintf("Serie che inizia prima del 1951 %s", ID)
        ## questo nn verrà mai eseguito
        x <- subset(x, x >= YYYY)
    } else {
        return()
    }

    x.temp <- x
    ## questo ci serve solo per il grafico (vogliamo mantenere gli NA)

    index <- which(is.na(x$prcp))
    if (length(index)) 
        x <- x[-index, ]

    yy.index <- matrix(1:length(x$prcp), ncol = 1)
    x <- cbind(yy.index, x)


    ## theil sen
    x.zyp <- zyp.sen(prcp ~ yy, dataframe = x)
    zyp.coef <- coefficients(x.zyp)

    zyp.int <- confint.zyp(x.zyp)

    plot(x.temp$yy, x.temp$prcp, type = "l", xlab = "Years", ylab = "mm", main = paste0("Station -- ", 
        ID))
    points(x.temp$yy, x.temp$prcp, pch = 21, bg = "red")
    abline(a = zyp.coef["Intercept"], b = zyp.coef["yy"])

    ## trend in the location parameter
    annual.fit.trend <- tryCatch({
        fevd(x = prcp, data = x, type = "GEV", location.fun = ~yy.index, use.phi = F, 
            units = "mm", time.units = "year")
    }, error = function(err) {
        print(paste(ID, "gev", sep = " "))
        print(err)
        my.logRatio <- c(ID = ID, DEV = NA, CHI = NA, MU = NA)
        return(my.logRatio)
    })


    parametri.trend <- round(annual.fit.trend$results$par[c(2, 3, 4)], 2)

    annual.fit.notrend <- 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)
        my.logRatio <- c(ID = ID, DEV = NA, CHI = NA, MU = NA)
        return(my.logRatio)
    })

    parametri.notrend <- round(annual.fit.notrend$results$par[c(1, 2, 3)], 2)

    print("Annual fit with trend in location parameter")
    plot(annual.fit.trend, , rperiods = tempi.ritorno)

    print("Annual fit with no trend in location parameter")
    plot(annual.fit.notrend, rperiods = tempi.ritorno)


    ## quale modello funziona meglio? loglikelihood ratio test
    lr.ris <- lr.test(annual.fit.notrend, annual.fit.trend, alpha = 0.05, df = 1)

    my.logRatio <- c(ID = ID, DEV = round(lr.ris$statistic, 2), PVAL = lr.ris$p.value, 
        ZYP.SLOPE = zyp.coef["yy"], ZYP.LW = zyp.int[2], ZYP.UP = zyp.int[4])


    list(my.logRatio)


}  ##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)

## parametri
oldPar <- par()

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

## Stazione: 6215

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6220

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6224

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6227

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6234

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6247

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6283

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6293

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6295

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6300

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6301

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6307

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6309

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6313

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6250

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6260

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6242

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6243

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6268

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6272

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6278

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6206

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6225

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1

## Stazione: 6261

plot of chunk unnamed-chunk-1

## [1] "Annual fit with trend in location parameter"

plot of chunk unnamed-chunk-1

## [1] "Annual fit with no trend in location parameter"

plot of chunk unnamed-chunk-1


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

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