Отгадка, метод EVT GPD Peak Over Thresholds действительно не работает. Даже при сэмпле 100к StudentT(df=4) он в половине случаев не может точно определить экспоненту хвоста, у него системная ошибка, биас, MLE фиттинг систематически недооценивает хвост, и получается как на картинке (это лог лог CDF правого хвоста, >0.99 квантили).
https://imgur.com/a/enrinVOФиттинг сделан верно, получилась некая "средняя" (точнее 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)
04.08.2025, 09:10 
(Hansen variant) позволяла бы разные степени для левого и правого 
