In MartinRoth/regthresh: Regional Threshold Analysis

library(regthresh)
library(copula)
library(rgdal)
library(extRemes)
library(potreg)
library(ggmap)
library(data.table)
library(doParallel)
registerDoParallel()

sites <- 16
years <- 50
blockLength <- 92
days  <- years * blockLength
probs <- seq(.9, .995, by = 0.001)
retPeriods <- c(5, 50, 100)

xpar <- list()
xpar$S <- sites
xpar$T <- days
xpar$fpar <- function(p, xpar) {
  loc   <- xpar$loc
  scale <- p[1] * loc
  shape <- matrix(p[2], xpar$T, xpar$S)
  list(loc = loc, scale = scale, shape = shape)
}
xpar$nOfP <- 2
xpar$start <- c(0.5, 0)

tau     <- 0.95
epsilon <- 0.25 
gamma   <- 0.5
xi      <- 0.15
kappasV <- rbeta(xpar$S, 2, 5) + 0.5
kappasP <- rep(calculateCorrespondingWeibullShape(0.15, tau), xpar$S)
betas   <- runif(16, 2, 4)

mParamsV <- setMarginalParameters(tau = tau, epsilon = epsilon, betas = betas
    ,kappas = kappasV, gamma = gamma, xi = xi, xpar = xpar) 
mParamsP <- setMarginalParameters(tau = tau, epsilon = epsilon, betas = betas
    ,kappas = kappasP, gamma = gamma, xi = xi, xpar = xpar)

x <- seq(0, 50, by = 0.1)
plot(x, pDW(x, mParamsP[[1]]), type = "l")
for(i in 2 : 5) lines(x, pDW(x, mParamsV[[i]]), lty = i)
abline(h = 0.95, col = 2, lty = 3)

x <- seq(12, 16, by = 0.1)
plot(x, dDW(x, mParamsP[[1]]), type = "l", ylim = c(0, 0.025))
lines(x, dDW(x, mParamsV[[1]]), col = 2)
for(i in 2 : 5) lines(x, dDW(x, mParamsV[[i]]), lty = i)

x <- seq(9, 20, by = 0.1)
plot(x, hazDW(x, mParamsP[[1]]), type = "l", ylim = c(0, .3))
lines(x, hazDW(x, mParamsV[[1]]), col = 2)
for(i in 2 : 5) lines(x, hazDW(x, mParamsV[[i]]), lty = i)

copulaInd <- normalCopula(0, sites)
dataVInd  <- generateDWData(mParamsV, copulaInd, xpar)
dataPInd  <- generateDWData(mParamsP, copulaInd, xpar)

Plot area of Vallei en Veluwe waterboard.

waterschap.rg      <- readOGR("./inst/extData/Waterschapsgrenzen_28992/","AU_AdministrativeUnitPolygon")
waterschap.sp      <- spTransform(waterschap.rg, CRS("+proj=longlat"))
ValleiEnVeluwe.sp  <- subset(waterschap.sp, waterschap.sp$name == "Vallei & Veluwe")
ValleiEnVeluwe.dat <- fortify(ValleiEnVeluwe.sp) 

basemap <- ggmap(get_map(location = c(lon = 5, lat = 52.5),
                    color = "color",
                    source = "stamen",
                    maptype = "watercolor",
                    zoom = 7))

focusBasemap <- ggmap(get_map(location = c(lon = 5.75, lat = 52.25),
                     color = "color",
                     source = "stamen",
                     maptype = "watercolor",
                     zoom = 9))

Load rain data

prDat <- knmiR::HomogenPrecip(ValleiEnVeluwe.sp, period="1951/2015")
prDat[, year := year(date)]
prDat[, month := month(date)]
prDat <- prDat[month %in% c(6,7,8)]

overview <- basemap + geom_path(aes(x = long, y = lat), col = "red", size = 1, data = ValleiEnVeluwe.dat)
print(overview)
setkey(prDat, stationId)
focusview <- focusBasemap + 
  geom_path(aes(x = long, y = lat), col = "red", size = 1, data = ValleiEnVeluwe.dat) +
  geom_point(aes(x = lon, y = lat), col = "black", shape = 3,
             data = knmiR::stationMetaData[.(prDat[, unique(stationId)]), ]) +
  geom_point(aes(x = lon, y = lat), col = "blue", shape = 4,
             data = knmiR::stationMetaData[stationName=="Putten", ])
print(focusview)

summary(prDat)

prDat[, jdays := julian(date, origin=as.Date("1951-01-01"))]
setkey(prDat, stationId)

#prDat[, prDcl2 := NA_real_]
#system.time(prDat[, prDcl := potreg::declusterData(jdays, pr, 1), by = stationId])
#system.time({
prDcl <- foreach(stnId=unique(prDat[["stationId"]]), .combine = "c") %dopar% {
  prDat[.(stnId), declusterData(jdays, pr, 1)]
}
prDat[, prDcl := prDcl]
#})

ObtainTSPlot <- function(x, tau) {
  tmp <- foreach(q = iter(tau), .combine = "rbind") %dopar% {
    data.table(tau = q, shape = fevd(x, threshold = quantile(x, q), type="GP")$results$par[2])
  }
  return(tmp)
}

setQuantiles <- function(dat, tau) {
  return(matrix(apply(dat, 2, quantile, tau),
                nrow=nrow(dat), ncol=ncol(dat), byrow=TRUE))
}

getExceedances <- function(data, u) {
  index        <- data > u
  data[!index] <- NA
  return(data)
}

getExcesses <- function(data, u) {
  excesses <- getExceedances(data - u, 0)
  return(excesses)
}

mTS1 <- prDat[stationId == knmiR::stationMetaData[stationName=="Putten", stationId], ObtainTSPlot(prDcl, probs)]
#plot(probs[-(93:96)], mTS1[-(93:96)], xlab = "", ylab = "")
#ggplot(mTS1, aes(x = tau, y = shape)) + geom_line() + xlim(c(0.9, 0.985)) + ylim(c(0, 0.3))

mSIVInd <- prDat[, ObtainTSPlot(prDcl, probs), by = stationId]
mSIVInd <- mSIVInd[, .(type = "individual", shape = mean(shape)), by = tau]
#ggplot(mSIVInd, aes(x = tau, y = shape)) + geom_line() + xlim(c(0.9, 0.985)) + ylim(c(0.025, 0.2))

tmpData <- as.matrix(dcast(prDat, date ~ stationId, value.var = 'prDcl')[, -1, with=FALSE])
tmpXpar <- list()
tmpXpar$S <- ncol(tmpData)
tmpXpar$T <- nrow(tmpData)
fpar <- function(p, xpar) {
  loc = xpar$loc
  scale = p[1]*loc
  shape = matrix(p[2], nrow = tmpXpar$T, ncol = tmpXpar$S)
  list(loc = loc, scale = scale, shape = shape)
}

mTSCV <- foreach(q = iter(probs), .combine = "rbind") %dopar% {
  tmpXpar$loc <- setQuantiles(tmpData, q)
  excess <- getExcesses(tmpData, tmpXpar$loc)
  tmp <- fitGpdFlex(excess, tmpXpar, fpar, hessian = TRUE,
                numberOfParameters = 2, start = c(0.5, 0),  method = "BFGS",
                control = list(maxit = 5000, ndeps = c(1e-03, 1e-05)))
  return(data.table(type = "joint", tau = q, shape = tmp$estimate[2]))
}

ggplot(rbind(mSIVInd, mTSCV), aes(x = tau, y = shape, col = type)) + 
  geom_line() + xlim(c(0.9, 0.985)) + ylim(c(0.025, 0.2))

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

When you click the Knit button a document will be generated that includes both content as well as the output of any embedded R code chunks within the document. You can embed an R code chunk like this:

MartinRoth/regthresh documentation built on May 7, 2019, 3:40 p.m.