Chebyshev, Extreme Value Theory ve Peak Over Threshold: Uç Olayların İstatistiği

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import genpareto, kstest
 
def mean_excess_plot(x, thresholds=None):
    """
    Mean excess plot — eşik seçimi için.
    """
    x = np.asarray(x)
    if thresholds is None:
        thresholds = np.linspace(np.percentile(x, 70), np.percentile(x, 99), 50)
    means = []
    for u in thresholds:
        excess = x[x > u] - u
        if len(excess) < 5:
            means.append(np.nan)
        else:
            means.append(excess.mean())
    return thresholds, np.array(means)
 
def fit_pot(x, threshold):
    """
    Generalized Pareto Distribution fit (above threshold).
    """
    excess = x[x > threshold] - threshold
    if len(excess) < 10:
        raise ValueError("Çok az excess gözlem")
    # scipy.stats.genpareto: (c, loc, scale) — c = ξ
    shape, loc, scale = genpareto.fit(excess, floc=0)
    return {'shape': shape, 'scale': scale, 'n_excess': len(excess)}
 
def pot_var(x, threshold, p=0.999, gpd_params=None):
    """
    POT ile VaR hesaplaması.
    p = güvence seviyesi (örn. 0.999 = %99.9)
    """
    if gpd_params is None:
        gpd_params = fit_pot(x, threshold)
    n = len(x)
    Nu = gpd_params['n_excess']
    xi = gpd_params['shape']
    sigma = gpd_params['scale']
 
    if abs(xi) < 1e-6:
        var = threshold - sigma * np.log((n / Nu) * (1 - p))
    else:
        var = threshold + (sigma / xi) * (((n / Nu) * (1 - p)) ** (-xi) - 1)
    return var
 
# Örnek: heavy-tail simülasyon
np.random.seed(42)
normal_part = np.random.normal(50, 5, 9000)
heavy_tail = np.random.pareto(a=2.5, size=1000) * 20 + 50  # heavy tail
data = np.concatenate([normal_part, heavy_tail])
 
# Mean excess plot
ths, mes = mean_excess_plot(data, thresholds=np.linspace(60, 120, 50))
plt.plot(ths, mes, 'o-')
plt.xlabel('Threshold u')
plt.ylabel('Mean Excess')
plt.title('Mean Excess Plot')
plt.grid()
plt.savefig('reports/mean_excess.png', dpi=120)
 
# u=80'i seç (görsele bakarak)
u = 80
gpd = fit_pot(data, threshold=u)
print(f"GPD fit: shape={gpd['shape']:.3f}, scale={gpd['scale']:.3f}, N_u={gpd['n_excess']}")
 
# 99.9% VaR
var_999 = pot_var(data, threshold=u, p=0.999)
print(f"VaR (99.9%): {var_999:.2f}")
# Bu eşiği aşan gözlemler 'anomali' kabul edilir.

# SPOT pseudo-code (Siffer et al., 2017)
class SPOT:
    """
    Streaming Peak Over Threshold — online anomaly detection.
    """
    def __init__(self, q=0.001, n_init=1000):
        self.q = q              # anomaly olasılığı (örn. 0.001)
        self.n_init = n_init    # warm-up veri sayısı
        self.threshold = None   # mevcut eşik (zT)
        self.peaks = []         # eşiği aşan excess'ler
        self.t = None           # GPD eşiği (u)
 
    def initialize(self, x_init):
        """Warm-up: ilk n_init veriyle GPD'yi fit et."""
        self.t = np.percentile(x_init, 98)
        excess = x_init[x_init > self.t] - self.t
        self.shape, _, self.scale = genpareto.fit(excess, floc=0)
        self.peaks = list(excess)
 
        # Anomaly threshold
        n = len(x_init)
        Nt = len(excess)
        self.threshold = self._compute_threshold(n, Nt)
 
    def _compute_threshold(self, n, Nt):
        xi, sigma = self.shape, self.scale
        if abs(xi) < 1e-6:
            return self.t - sigma * np.log((n / Nt) * self.q)
        return self.t + (sigma / xi) * (((n / Nt) * self.q) ** (-xi) - 1)
 
    def update(self, x_new):
        """Bir yeni gözlem geldiğinde güncelle."""
        is_anomaly = x_new > self.threshold
        if x_new > self.t and not is_anomaly:
            # Eşik üstü ama anomaly değil — GPD'yi güncelle
            self.peaks.append(x_new - self.t)
            self.shape, _, self.scale = genpareto.fit(self.peaks, floc=0)
            n = len(self.peaks) * 50  # estimated n
            self.threshold = self._compute_threshold(n, len(self.peaks))
        return is_anomaly
 
# Kullanım (kavramsal — production'da daha karmaşık)
spot = SPOT(q=0.001, n_init=1000)
spot.initialize(data[:1000])
 
for x in data[1000:]:
    if spot.update(x):
        print(f"Anomaly: {x}")