Estimating Heavy-Tail Exponents with Hill and Pareto Methods

# EVT - How to estimate tail exponent with good precision? Hill, GPT EVT, LeastSquares


# EVT - How to estimate tail exponent with good precision? Hill, GPT EVT, LeastSquares












I generated 30 samples of 20k points of StudentT(df=4) and tried various tail estimators for df. All results are terrible.

Line Color - type of estimator, x - the estimator treshold value, each line - separate simulation (sample). Expected correct result - constant line with y = 4.

I assumed the GPD (Generalised Pareto Distribution estimator from EVT) would be the best one, yet it seems to be the worst (blue lines), or maybe I implemented it wrongly.

The code for GPD estimator:

```
def student_sample(df, n, seed=None):
  rng = np.random.default_rng(seed)
  return rng.standard_t(df, size=n)

def estimate_gpd_mle(x, threshold, init=(1/5.0, 1.0), min_exceed=5):
  # Generalized Pareto Distribution (GPD) MLE estimation
  y = x[x > threshold] - threshold
  if y.size < min_exceed: return None

  def nll(params):
    ξ, β = params
    if β <= 0: return np.inf
    base = 1 + ξ * y / β
    if np.any(base <= 0): return np.inf
    pdf = (1 / β) * base ** -(1/ξ + 1)
    if np.any(pdf <= 0) or np.any(~np.isfinite(pdf)): return np.inf
    return -np.sum(np.log(pdf))
  res = minimize(nll, x0=init, bounds=((1/20, 1/2), (1e-3, 10)))
  if not res.success: return None
  ξ, β = res.x
  return (1/ξ, β)  # return df = 1/xi
```

The other approaches Hill (black), Least Squares (green) - seems also work poorly, failing to find correct df=4 and producing biased estimation ~3.5. The full distribution MLE fit (red) produce the best result, but it's not a tail estimator.

There's advice to reject top 5-15 points, but it doesn't change the picture much (the chart below). Also, for Hill estimator (black lines) the selection of threshold is "when the chart starts to stabilise" - good look finding that - it seems to be quite arbitrary, you can choose any point marked with yellow circle, they all could be considered in some sense as "starting to stabilise" (especially if red lines are removed and we don't know that true df = 4).

So, is there a way to more or less reliably estimate the tails? I would like to find out if return distributions for T=30,180,365 have different tail exponents. And I can't use full Distribution MLE fit, because StudentT or SkewStudentT or GenHyperbolic all distort the tails and also tail affected by the skew, so the full distribution MLE fit won't produce the true df tail parameter.

I guess it's possible to use Hill estimator, if tail exponent for left and right tails are different - it should show the difference (even if it won't allow to get the absolute value of df with reasonable precision).

UPDATE

A log log plot of right tail for 6 samples, with true lines df=3.5,4,4.5. Seems like the visual approach, should produce no worse result than the Hill

Full code to reproduce the chart

```
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t
from scipy.optimize import minimize_scalar, minimize

def student_sample(df, n, seed=None):
  rng = np.random.default_rng(seed)
  return rng.standard_t(df, size=n)

def estimate_hill(x):
  x = np.sort(x)[::-1]
  logx = np.log(x)
  sumk = np.cumsum(logx)
  k = np.arange(1, len(x)+1)
  denom = (sumk/k) - logx
  return np.array([1/d if d != 0 else np.nan for d in denom])

def fit_mle_student(x):
  def nll(df):
    return np.inf if df <= 2 else -np.sum(t.logpdf(x, df, loc=0, scale=1))
  res = minimize_scalar(nll, bounds=(2.01, 100), method='bounded')
  if not res.success:
    raise RuntimeError("MLE fit failed")
  return res.x

def estimate_gpd_mle(x, threshold, init=(1/5.0, 1.0), min_exceed=5):
  # Generalized Pareto Distribution (GPD) MLE estimation
  y = x[x > threshold] - threshold
  if y.size < min_exceed: return None

  def nll(params):
    ξ, β = params
    if β <= 0: return np.inf
    base = 1 + ξ * y / β
    if np.any(base <= 0): return np.inf
    pdf = (1 / β) * base ** -(1/ξ + 1)
    if np.any(pdf <= 0) or np.any(~np.isfinite(pdf)): return np.inf
    return -np.sum(np.log(pdf))
  res = minimize(nll, x0=init, bounds=((1/20, 1/2), (1e-3, 10)))
  if not res.success: return None
  ξ, β = res.x
  return (1/ξ, β)  # return df = 1/xi

def estimate_gpd_mle_and_prepare_data_for_plot(x, K, count=20, init=(0.2, 1.0)):
  # Ignore, nothing interesting, just some data preparation for the plot
  x_sorted = np.sort(x)
  ranks = np.linspace(10, K, count).astype(int)
  df_estimates = np.full_like(ranks, np.nan, dtype=float)
  current_init = init

  for i, k in enumerate(ranks):
    threshold = x_sorted[-k]
    res = estimate_gpd_mle(x, threshold, init=current_init)
    if res is None: continue
    df, beta = res
    df_estimates[i] = df
    # warm start next with xi,beta
    current_init = (1/df, beta)

  return ranks, df_estimates

def estimate_tail_ls(x):
  # Least Squares estimation of the tail index
  n = len(x)
  if n < 3:
    return None
  x_sorted = np.sort(x)
  logx = np.log(x_sorted)
  logp = np.log(np.arange(n, 0, -1)) - np.log(n + 1)
  A = np.vstack([logx, np.ones_like(logx)]).T
  try:
    slope, _ = np.linalg.lstsq(A, logp, rcond=None)[0]
    alpha = -slope
    return alpha if alpha > 0 else None
  except:
    return None

def estimate_tail_ls_and_prepare_data_for_plot(x, K, count=20):
  # Ignore, nothing interesting, just some data preparation for the plot
  x_sorted = np.sort(x)
  ranks = np.linspace(10, K, count).astype(int)
  df_estimates = np.full_like(ranks, np.nan, dtype=float)

  for i, k in enumerate(ranks):
    tail = x_sorted[-k:]
    df = estimate_tail_ls(tail)
    if df is not None:
      df_estimates[i] = df

  return ranks, df_estimates

# # Parameters
N = 30
n_obs = 20_000
q = 0.97

plt.figure(figsize=(8, 6))
for _ in range(N):
  x = student_sample(4, n_obs)

  # Hill estimates
  k = int((1 - q) * len(x))
  tail = np.sort(x)[-k:][::-1]
  hill = estimate_hill(tail)
  plt.plot(hill, linewidth=1, alpha=0.5, color='black')

  # Exact t MLE
  mle_df = fit_mle_student(x)
  plt.axhline(mle_df, linewidth=1, alpha=0.5, color='red')

  # GPD estimates across multiple order-statistic thresholds
  ranks, dfs = estimate_gpd_mle_and_prepare_data_for_plot(x, K=k, count=20)
  plt.plot(ranks, dfs, marker='o', linestyle='-', markersize=3, alpha=0.5, color='blue', linewidth=1, )

  # LS estimates across multiple order-statistic thresholds
  ranks, dfs = estimate_tail_ls_and_prepare_data_for_plot(x, K=k, count=20)
  plt.plot(ranks, dfs, marker='o', linestyle='-', markersize=3, alpha=0.5, color='green', linewidth=1)

plt.ylim(2, 6)
plt.xlabel('Order-stats rank k')
plt.ylabel('Tail index estimate')
plt.title(f'Hill (black), MLE (red), GPD (blue), LS (green) over {N} samples. True ν = 4')
plt.show()
```

UPDATE2

R code from POT library, same terrible result.

R code

```
library(POT)
set.seed(1)

estimate_gpd_mle <- function(x, k) {
  u <- sort(x, decreasing = TRUE)[k]
  xi <- coef(fitgpd(x, threshold = u, est = "mle"))["shape"]
  as.numeric(1 / xi)
}

estimate_gpd_mle_for_plot <- function(x, tail_q, count = 20) {
  K <- floor(tail_q * length(x))
  ranks <- as.integer(seq(10, K, length.out = count))
  df_estimates <- rep(NA_real_, length(ranks))
  for (i in seq_along(ranks)) {
    df_estimates[i] <- tryCatch(
      estimate_gpd_mle(x, ranks[i]),
      error = function(e) NA_real_
    )
  }
  data.frame(k = ranks, df = df_estimates)
}

# estimate_gpd_mle_for_plot(x, K = 200, count = 5)

plot_gpd_estimates <- function(n_points, n_samples, tail_q, drop_top_n, true_nu) {
  xlim_max <- floor(tail_q * n_points)
  plot(NULL, xlim = c(10, xlim_max), ylim = c(2, 8),
       xlab = "Order-stats rank k", ylab = "Tail index estimate",
       main = sprintf("Hill (black), MLE (red), GPD (blue), LS (green) over %d samples. True \u03BD = %d",
                      n_samples, true_nu))

  for (i in seq_len(n_samples)) {
    x <- rt(n_points, df = true_nu)
    if (drop_top_n > 0) {
      x <- sort(x, decreasing = TRUE)
      if (length(x) > drop_top_n) x <- x[(drop_top_n + 1):length(x)]
    }

    df_plot <- estimate_gpd_mle_for_plot(x, tail_q = tail_q, count = 20)
    lines(df_plot$k, df_plot$df, type = "b", pch = 16, cex = 0.6, lwd = 1, col = "blue", lty = 1)
  }

  abline(h = true_nu, col = "red", lty = 2)
  legend("topright",
         legend = c("GPD MLE runs", sprintf("True \u03BD = %d", true_nu)),
         col = c("blue", "red"), lty = c(1, 2), pch = c(16, NA), bty = "n")
}

plot_gpd_estimates(n_points = 20000, n_samples = 3, tail_q = 0.01, drop_top_n = 0, true_nu = 4)
```

## Answer by Con Fluentsy (score 0)

https://quant.stackexchange.com/a/84152

It is not suprising you endedup with a mess, you mixed metaphors as they say in english teaching, stick to one framework at a time, do not mix methods. Here is an example of a simple POT model working fine:

```
library(evir)
library(quantmod)

# ---- Get Financial Data ----
ticker <- "QQQ"
beg <- "2000-01-01"
ed <- "2025-03-05"

# Get stock data
getSymbols(ticker, from=as.Date(beg), to=as.Date(ed))
price.mat <- Cl(QQQ)

# Compute daily log returns
QQQ_r <- dailyReturn(price.mat, type = "log")

# ---- Fit Generalized Pareto Distribution (GPD) ----
fit <- gpd(QQQ_r, threshold = quantile(QQQ_r, 0.90))  # Using 90th percentile threshold

# ---- Generate Tail Plot Object ----
tp <- tailplot(fit)  

# ---- Estimate Extreme Quantiles ----
# Estimate the 99.9% quantile with 95% confidence interval
extreme_quantile <- gpd.q(tp, pp = 0.999, ci.p = 0.95)  
print(extreme_quantile)

# Compare with direct sample quantile estimation
sample_quantile <- quantile(QQQ_r, probs = 0.999)
print(sample_quantile)

# ---- Compute Expected Shortfall (CVaR) Beyond the 99% Quantile ----
expected_shortfall <- gpd.sfall(tp, 0.99)  
print(expected_shortfall)

last(QQQ$QQQ.Close)
last_trade <- tail(QQQ$QQQ.Close, 1)
expected_shortfall * c(last_trade,last_trade, last_trade)

hist(QQQ_r, breaks = 200, xlim = c(min(QQQ_r), max(QQQ_r)))
abline(v = 0, col = "red", lwd = 2)
abline(v = quantile(QQQ_r, 0.10), col = "blue", lwd = 2, lty = 2)
  legend("topright", legend = c("Zero Line", "10th Percentile"), 
       col = c("red", "blue"), lwd = 2, lty = c(1, 2))
```
```