Отгадка, EVT GPD POT оценка хвостов действительно не работает

Я столкнулся с проблемой — самый точный (в теории) метод оценки экспоненты хвоста распределения давал совершенно дикие результаты smart-lab.ru/blog/1188320.php Мне кажется я нашел причину почему. Он действительно ошибается.

Даже при сэмпле 100к StudentT(df=4) он в половине случаев не может точно определить экспоненту хвоста, у него сильная системная ошибка, биас, MLE фиттинг систематически недооценивает хвост, и получается как на картинке (это лог лог CDF правого хвоста, >0.99 квантили). 


Фиттинг сделан верно, получилась некая «средняя» (точнее MLE) линия проходящая по всем точкам, и которая систематически недооценивает экспоненту.

И становится понятно эмпирическое правило «выкинуть последние 5-10 точек» — точнее, это правило неверное, если просто выкинуть случайно несколько последних точек результат буде вообще дикий, видимо в этом правиле есть упущенное дополнение, что выкидывать нужно не просто случайно точки, а только последние точки которые пошли вниз.

P.S. Как это сделать автоматически — может выкинуть последние точки, как вариант — нужно итерационно делать фиттинг, и выкидывать последнюю точку до тех пор пока df растет, а как только он упал сразу остановиться. Или, попробовать другие эстиматоры кроме MLE…

R если кому интересно посмотреть

library(POT)
set.seed(1)

estimate_gpd_mle <- function(x, u) {
  coef_vec <- coef(fitgpd(x, threshold = u, est = "mle"))
  list(
    scale = as.numeric(coef_vec["scale"]),
    shape = as.numeric(coef_vec["shape"])
  )
}

estimate_gpd_mle_df <- function(x, u) {
  as.numeric(1 / estimate_gpd_mle_df(x, u))
}

x <- rt(100000, df = 4)
quantiles <- seq(0.95, 0.995, length.out = 50)
dfs <- sapply(quantiles, function(q) {
  u <- quantile(x, q)
  estimate_gpd_mle(x, u)
})


plot(quantiles, dfs, type = "b", pch = 19, col = "blue",
     xlab = "Quantile threshold (q)",
     ylab = "Estimated degrees of freedom (1/xi)",
     main = "Tail Index Estimation via GPD", log = "y", ylim = c(2, 8))
abline(h = 4, col = "red", lty = 2)  # true df
legend("topright", legend = "True df = 4", col = "red", lty = 2, bty = "n")


plot_tail_ll <- function(x, q, true_df, fit) {
  # Extract threshold
  u <- quantile(x, q)
  
  # Empirical tail
  tail_x <- x[x > u]
  tail_x_sorted <- sort(tail_x)
  
  n <- length(tail_x_sorted)
  ccdf <- rev(seq_along(tail_x_sorted)) / n
  
  # Plot empirical CCDF
  plot(tail_x_sorted, ccdf, log = "xy", type = "p", pch = 19, cex = 0.5,
       xlab = "x (log scale)",
       ylab = "P(X > x | X > u) (log scale)",
       main = paste("Log-Log Tail CCDF (q =", q, ")"))
  
  # Fitted GPD CCDF
  xi <- fit$shape
  sigma <- fit$scale
  
  x_grid <- seq(min(tail_x_sorted), max(tail_x_sorted), length.out = 200)
  excesses <- x_grid - u
  ccdf_gpd <- (1 + xi * excesses / sigma) ^ (-1 / xi)
  
  # Add fitted GPD line
  lines(x_grid, ccdf_gpd, col = "red", lwd = 2)
  
  # Add true Student-t tail
  ccdf_true <- 1 - pt(x_grid, df = true_df)
  lines(x_grid, ccdf_true / mean(x > u), col = "black", lty = 2, lwd = 1)  # conditional CCDF
  
  # Annotate
  legend("topright",
    legend = c(
      sprintf("GPD fit: df = %.4f, scale = %.4f", 1/xi, sigma),
      sprintf("True Student-t df = %d", true_df)
    ), col = c("red", "black"), lty = c(1, 2), bty = "n")
}


fit <- estimate_gpd_mle(x, quantile(x, 0.99))
plot_tail_ll(x, 0.995, 4, fit)

# x_trunc <- sort(x)[1:(length(x) - 5)]
# fit <- estimate_gpd_mle(x_trunc, quantile(x, 0.99))
# plot_tail_ll(x_trunc, 0.995, 4, fit)