Расчет уровней возврата на основе GPD в разных R-пакетах

Question

Я выполняю анализ экстремальных значений для метеорологических данных, а точнее для данных о осадках, доступных в мм/д. Я использую пороговый избыточный подход для оценки параметров обобщенного распределения Парето с методом максимального правдоподобия.

Цель состоит в том, чтобы рассчитать несколько уровней возврата (т. Е. События 2, 5, 10, 20, 50, 100 год) для суточного осаждения.

Хотя код R работает нормально, мне интересно, почему я получаю явно разные результаты при расчете уровней возврата на основе квантилей установленного GPD с различными пакетами. Несмотря на то, что оценочные параметры GPD почти одинаковы в каждой упаковке, квантилисты сильно различаются.

Пакеты, которые я использовал: ismev, extRemes, evir и POT.

Я полагаю, что разные оценки параметров GPD связаны с различными процедурами расчета, но я не понимаю, почему вычисление квантилей сильно различается в зависимости от разных пакетов.

В то время как lmom, evir и POT возвращают одинаковые значения, возвращаемый период, полученный из пакета extRemes, отличается от других результатов.

# packages 
library(ismev) 
library(extRemes) 
library(evir) 
library(POT) 
library(lmom) 

th <- 50 

# sample data: 
potvalues <- c(
    58.5,44.2,49.6,59.3,48.3,60.9,94.5,47.1,45.3,57.6,48.2,46.2,44.2,50.6,42.1,52.7,80.9, 
    58.5,51.3,48.4,51.7,71.9,60.1,64.4,43.5,55.5,49.3,58.2,47.5,43.7,45.2,52.8,42.2,46.4, 
    96.1,47.5,50.1,42.4,60.9,72.6,51.6,59.4,80.5,63.7,59.9,45.0,66.7,47.6,53.3,43.1,51.0, 
    46.2,53.6,59.8,51.7,46.7,42.6,44.5,45.0,50.0,44.0,89.9,44.2,47.8,53.3,43.0,55.7,44.6, 
    44.6,54.9,45.1,43.9,78.7,45.5,64.0,42.7,47.4,57.0,105.4,64.3,43.2,50.4,80.2,49.9,71.6, 
    47.4,44.1,47.6,55.2,44.4,78.6,50.8,42.4,47.1,43.5,51.4) 

#------------------------------------------------------------------------------------------# 

# MLE Fitting of GPD - package extRemes 

# fit gpd 
pot.ext <- fevd(potvalues, method = "MLE", type="GP", threshold=th) 

# return levels: 
rl.extremes <- return.level(pot.ext, conf = 0.05, 
          return.period= c(2,5,10,20,50,100)) 
rl.extremes <- as.numeric(rl.extremes) 

#------------------------------------------------------------------------------------------# 

# MLE Fitting of GPD - package ismev 

pot.gpd <- gpd.fit(potvalues, threshold=th) 

s1 <- quagpa(f=.99, para=c(pot.gpd$threshold, pot.gpd$mle[1],-pot.gpd$mle[2])) # 100 
s2 <- quagpa(f=.98, para=c(pot.gpd$threshold, pot.gpd$mle[1],-pot.gpd$mle[2])) # 50 
s3 <- quagpa(f=.95, para=c(pot.gpd$threshold, pot.gpd$mle[1],-pot.gpd$mle[2])) # 20 
s4 <- quagpa(f=.90, para=c(pot.gpd$threshold, pot.gpd$mle[1],-pot.gpd$mle[2])) # 10 
s5 <- quagpa(f=.80, para=c(pot.gpd$threshold, pot.gpd$mle[1],-pot.gpd$mle[2])) # 5 
s6 <- quagpa(f=.50, para=c(pot.gpd$threshold, pot.gpd$mle[1],-pot.gpd$mle[2])) # 2 

rl.ismev <- c(s6, s5, s4, s3, s2, s1) 

#------------------------------------------------------------------------------------------# 

# MLE Fitting of GPD - package evir 

# fit gpd 
gpd.evir <- gpd(potvalues, threshold=th) 

# plot 
evirplot <- plot(gpd.evir) 
1 # Excess Distribution 
0 # exit 

x100 <- gpd.q(pp=.99, x=evirplot) # 100 
x050 <- gpd.q(pp=.98, x=evirplot) # 50 
x020 <- gpd.q(pp=.95, x=evirplot) # 20 
x010 <- gpd.q(pp=.90, x=evirplot) # 10 
x005 <- gpd.q(pp=.80, x=evirplot) # 5 
x002 <- gpd.q(pp=.50, x=evirplot) # 2 

rl.evir <- t(rbind(x002,x005,x010,x020,x050,x100)) 
rl.evir <- as.numeric(rl.evir[2,]) 

#------------------------------------------------------------------------------------------# 

# MLE Fitting of GPD - package POT 

gpd.pot <- fitgpd(potvalues, threshold=th) 
quant = c(0.50, 0.80, 0.90, 0.95, 0.98, 0.99) 
rtp <- c(2,5,10,20,50,100) 

retvec <- vector() 
for (i in quant){ 
    x <- POT::qgpd(i, loc = th, scale = as.numeric(gpd.pot$param[1]), 
      shape = as.numeric(gpd.pot$param[2])) 
    retvec <- c(retvec,x) 
} 

rl.pot <- retvec 

#------------------------------------------------------------------------------------------# 
# comparison of results - return periods 
result <- cbind(rl.extremes,rl.ismev, rl.evir, rl.pot) 
round(result, 2) 

#------------------------------------------------------------------------------------------# 
# comparison of estimated parameters 
param.extremes <- pot.ext$results$par # extremes 
param.ismev <- pot.gpd$mle # ismev 
param.evir <- c(gpd.evir$par.ests[2],gpd.evir$par.ests[1]) # evir 
param.pot <- gpd.pot$param # POT 

parameters <- cbind(param.extremes, param.ismev , param.evir, param.pot) 
round(parameters, 4) 

#------------------------------------------------------------------------------------------# 

Answer 1

Решение этой проблемы описано, например, в книге Coles (Введение в статистическое моделирование экстремальных значений, глава 4.3.3). Хотя уровни возврата для ГЭВ могут быть получены скорее непосредственно из его квантилей, так называемый коэффициент превышения (то есть количество событий в год или вероятность того, что событие превышает пороговое значение соответственно), необходимо учитывать при расчете уровней возврата для GP в рамках пика над пороговым приложением.

Уровень N-год возврата определяется

Таким образом, он не работает, чтобы получить значимые результаты для уровней возврата, когда просто вычисления квантилей для распределения ГП без учета скорости превышения. В пакете extRemes учитывается коэффициент превышения, в то время как значение по умолчанию для количества событий в год в пакетах POT и evir устанавливается равным 1, если не указано.

Answer 2

Различия могут также исходить из разных способов подгонки функции распределения к набору данных. У меня есть пакет на CRAN, сравнивающий GPD припадки (или, вернее, их квантильные оценки) несколько пакетов и методов R:

https://cran.r-project.org/web/packages/extremeStat/vignettes/extremeStat.html

Вы также можете использовать пакет для сравнения GPD с другими дистрибутивами.