# Extreme Value Theory (EVT) and Tail-Related Risk Measures

# 1. 尾部風險測度與核心目標

## 1.1 Value-at-Risk (VaR) 與 Expected Shortfall (ES) 定義

- **VaR (Value-at-Risk)：**  
  在給定置信水平 \\(\alpha\\)（例如 \\(\alpha=99.5\%\\)）下，VaR 是預計在特定時間內損失不會超過的水平。數學上，設 \\(L\\) 為損失（通常取正值表示損失大小），則 VaR 定義滿足：
  \\[
  P(L \leq VaR_\alpha) = \alpha
  \\]
  換句話說，如果一個投資組合的 \\(99.5\%\\) VaR 為 5%，則表示有 \\(99.5\%\\) 的概率，該損失不會超過 5%。

- **ES (Expected Shortfall) / CVaR (Conditional VaR)：**  
  ES 定義為在損失超過 VaR 的情況下的平均損失，即：
  \\[
  ES_\alpha = E[L \mid L > VaR_\alpha]
  \\]
  ES 描述了當損失超出 VaR 時，所預期的平均損失程度，能夠更全面地捕捉尾部風險，尤其在極端環境下更具意義。

## 1.2 核心目標

量化在極端情況下可能出現的損失，有助制定更保守的準備和風險控制策略。然而，傳統方法（如基於標準差的計算）常假設收益近似正態分佈，可能低估極端風險。為此，我們引入極值理論（EVT）來建模尾部行為，外推出極端損失的概率與幅度。

---

# 2. 極值理論（EVT）基本框架與主要方法

EVT 旨在描述隨機變量尾部的行為，不依賴整體分佈的假設。主要有兩種方法：

## 2.1 Block Maxima Method（區塊極值法）

- **基本思路：**  
  將數據 \\(X_1, X_2, \dots, X_N\\) 劃分為若干不重疊的區塊（例如每年或每月），並從每個區塊中提取極值 \\(M_i\\)（最大值或最小值）。

- **理論基礎：**  
  根據 [Fisher–Tippett–Gnedenko 定理](https://en.wikipedia.org/wiki/Fisher%E2%80%93Tippett%E2%80%93Gnedenko_theorem)，存在常數 \\(a_n>0\\) 和 \\(b_n\\)（分別為尺度與位置調整因子），使得當區塊長度 \\(n\\) 趨於無窮大時，
  
  \\[
  \frac{M_n - b_n}{a_n} \xrightarrow{d} G(z)
  \\]
  其中 \\(z = \frac{M_n - b_n}{a_n}\\) 為標準化後的極值，而 \\(G(z)\\) 的形式為[廣義極值分佈（GEV）](https://en.wikipedia.org/wiki/Generalized_extreme_value_distribution)的累積分佈函數 (CDF)：
  
  \\[
  G(z)=\exp\left\{-\left[1+\xi\left(\frac{z-\mu}{\sigma}\right)\right]^{-1/\xi}\right\},\quad 1+\xi\left(\frac{z-\mu}{\sigma}\right)>0
  \\]
  
  **符號說明：**
  - \\(M_n\\)：某區塊內的極值（例如最大損失）。
  - \\(a_n\\)：尺度調整因子（正標量），用於標準化極值的波動幅度。
  - \\(b_n\\)：位置調整因子，類似於均值調整。
  - \\(\mu\\)、\\(\sigma\\)、\\(\xi\\)：GEV 分佈的參數，其中 \\(\mu\\) 為位置參數、\\(\sigma>0\\) 為尺度參數、\\(\xi\\) 為形狀參數，描述尾部厚度。
  
- **應用注意事項：**  
  區塊極值法雖然理論上非常堅實，但每個區塊僅取一個極值，當樣本數有限時會浪費大量信息；而且區塊劃分（例如選擇每年、每月）的決策對結果有顯著影響。

## 2.2 Peaks Over Threshold (POT) Method（超過閾值法）

- **基本思路：**  
  選定一個較高閾值 \\(u\\)（例如數據的 95% 或 99% 分位數），並將所有超過該閾值的數據視為極端事件。定義超額數據：
  
  \\[
  Y = X - u,\quad \text{當 } X > u
  \\]

- **理論基礎：**  
  根據 [Pickands–Balkema–de Haan 定理](https://en.wikipedia.org/wiki/Pickands%E2%80%93Balkema%E2%80%93De_Haan_theorem)，當 \\(u\\) 足夠高時，超額數據的條件累積分佈
  
  \\[
  F_u(y)=\frac{F(u+y)-F(u)}{1-F(u)}
  \\]
  可近似為[廣義帕累托分佈（GPD）](https://en.wikipedia.org/wiki/Generalized_Pareto_distribution)的形式：
  
  \\[
  F_u(y) \approx 1-\left(1+\frac{\xi\,y}{\beta}\right)^{-1/\xi},\quad y\ge 0,\quad 1+\frac{\xi\,y}{\beta}>0
  \\]
  當 \\(\xi=0\\) 時退化為：
  
  \\[
  F_u(y) \approx 1-\exp\left(-\frac{y}{\beta}\right)
  \\]
  
  **符號說明：**
  - \\(u\\)：選定的高閾值。
  - \\(X\\)：原始數據（如損失）。
  - \\(Y = X - u\\)：超過閾值的部分，稱為超額數據。
  - \\(F(x)=P(X\le x)\\)：\\(X\\) 的累積分佈函數。
  - \\(\xi\\)：GPD 的形狀參數，描述尾部厚度（較大的 \\(\xi\\) 表示尾部較厚）。
  - \\(\beta\\)：GPD 的尺度參數，反映超額數據的變動規模。
  - \\(F_u(y)\\)：超額數據的條件累積分佈函數。

- **推導詳解：**

  1. **條件分佈定義：**  
     
     \\[
     F_u(y)=\frac{F(u+y)-F(u)}{1-F(u)}
     \\]
  
  2. **正則變化假設：**  
     假設 \\(X\\) 的尾部滿足[正則變化（Regular Variation；許多厚尾分佈都有的尾部特性，如：Pareto 或 Student's t）](https://en.wikipedia.org/wiki/Slowly_varying_function)，即：
     
     \\[
     1-F(x)=x^{-1/\xi}L(x)
     \\]
     其中 \\(L(x)\\) 是緩變函數，滿足對任何 \\(c>0\\)：
     
     \\[
     \lim_{x\to\infty}\frac{L(cx)}{L(x)}=1
     \\]
     這意味著當 \\(x\\) 趨於無窮大時，尾概率主要以冪律 \\(x^{-1/\xi}\\) 衰減。

  3. **推導條件分佈：**  
     當 \\(x=u+y\\) 時，
     
     \\[
     1-F(u+y) \approx (u+y)^{-1/\xi}L(u)
     \\]
     而
     
     \\[
     1-F(u) \approx u^{-1/\xi}L(u)
     \\]
     所以，
     
     \\[
     \frac{1-F(u+y)}{1-F(u)} \approx \left(\frac{u+y}{u}\right)^{-1/\xi} = \left(1+\frac{y}{u}\right)^{-1/\xi}
     \\]
     故：
     
     \\[
     F_u(y) \approx 1-\left(1+\frac{y}{u}\right)^{-1/\xi}
     \\]
     為了讓模型更靈活，引入尺度參數 \\(\beta\\)（通過擬合得到），最終得到：
     
     \\[
     F_u(y) \approx 1-\left(1+\frac{\xi y}{\beta}\right)^{-1/\xi}
     \\]
     當 \\(\xi=0\\) 時，則有：
     
     \\[
     F_u(y) \approx 1-\exp\left(-\frac{y}{\beta}\right)
     \\]
  
  4. **外推 VaR 與 ES：**  
     假設總樣本數為 \\(n\\)，超過閾值 \\(u\\) 的數據點數為 \\(n_u\\)，置信水平設為 \\(\alpha\\)（例如 99.5%），我們要求 \\(x_\alpha\\) 滿足：
     
     \\[
     P(X>x_\alpha) = 1-\alpha
     \\]
     設 \\(x_\alpha = u+y_\alpha\\)，則：
     
     \\[
     P(X > u+y_\alpha) \approx \frac{n_u}{n}\left(1+\frac{\xi y_\alpha}{\beta}\right)^{-1/\xi}
     \\]
     令上式等於 \\(1-\alpha\\)：
     
     \\[
     \frac{n_u}{n}\left(1+\frac{\xi y_\alpha}{\beta}\right)^{-1/\xi} = 1-\alpha
     \\]
     解得：
     
     \\[
     1+\frac{\xi y_\alpha}{\beta} = \left(\frac{n}{n_u}(1-\alpha)\right)^{-\xi}
     \\]
     進而，
     
     \\[
     y_\alpha = \frac{\beta}{\xi}\left[\left(\frac{n}{n_u}(1-\alpha)\right)^{-\xi}-1\right]
     \\]
     所以，
     
     \\[
     VaR_\alpha = x_\alpha = u+y_\alpha = u + \frac{\beta}{\xi}\left[\left(\frac{n}{n_u}(1-\alpha)\right)^{-\xi}-1\right]
     \\]
  
     ES 定義為在損失超過 VaR 的情況下的平均損失：
     
     \\[
     ES_\alpha = E[X\mid X>x_\alpha] = u + E[Y\mid Y>y_\alpha]
     \\]
     對於 GPD，經推導可得：
     
     \\[
     E[Y\mid Y>y_\alpha] = \frac{\beta + \xi y_\alpha}{1-\xi}
     \\]
     因此，
     
     \\[
     ES_\alpha = u + \frac{\beta}{1-\xi}\left[\left(\frac{n}{n_u}(1-\alpha)\right)^{-\xi}-1\right]
     \\]
     （該公式在 \\(\xi < 1\\) 時有效。）

---

# 3. 參數估計

## 3.1 Block Maxima Method
- 我們對每個區塊提取極值 \\(M_1, M_2, \dots, M_k\\)，然後利用最大似然估計（MLE）對這些數據直接擬合 GEV 分佈，估計出位置參數 \\(\mu\\)、尺度參數 \\(\sigma\\) 和形狀參數 \\(\xi\\)。
- 理論上存在標準化常數 \\(a_n\\) 和 \\(b_n\\) 使得 \\(z_i = \frac{M_i-b_n}{a_n}\\) 服從 GEV 分佈，但在 MLE 擬合中，我們直接用 \\(M_i\\) 進行擬合，這些常數效果被吸收到估計參數中。

## 3.2 POT 方法
- 對於超過閾值 \\(u\\) 的數據 \\(Y=X-u\\)，利用 MLE 擬合 GPD，得到形狀參數 \\(\xi\\) 與尺度參數 \\(\beta\\)。
- MLE 是最常用的方法，但也可以使用概率加權矩（PWM）。

---

# 4. Python 實作案例（基於 POT 方法）

以下代碼展示如何利用 POT 方法擬合 GPD 並計算 99.5% VaR 與 ES：

```python
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.stats import genpareto, t

# Step 1: Generate simulated daily returns with heavy tails (t distribution)
np.random.seed(42)
data = np.random.standard_t(df=3, size=1000) / 100.0  # returns ~ t(df=3) scaled by 1/100
returns = pd.Series(data)

# Step 2: Convert returns to losses (loss = -return)
# Here, positive values indicate losses.
losses = -returns

# Step 3: Set threshold as the 95th percentile of losses
threshold = np.percentile(losses, 95)
print("Threshold (95th percentile):", threshold)

# Step 4: Calculate exceedances (losses above threshold)
exceedances = losses[losses > threshold] - threshold
print("Number of exceedances:", len(exceedances))

# Step 5: Fit the Generalized Pareto Distribution (GPD) to exceedances using MLE
params = genpareto.fit(exceedances)
xi, loc, beta = params
print("Fitted GPD parameters:")
print("Shape (xi):", xi)
print("Scale (beta):", beta)

# Step 6: Calculate VaR and ES at 99.5% confidence level using POT formulas
alpha = 0.995
n = len(losses)          # total number of samples
n_exceed = len(exceedances)  # number of exceedances

# VaR formula:
var = threshold + (beta/xi) * (((n/float(n_exceed))*(1-alpha))**(-xi) - 1)
print("Estimated 99.5% VaR:", var)

# ES formula:
es = threshold + (beta/(1-xi)) * (((n/float(n_exceed))*(1-alpha))**(-xi) - 1)
print("Estimated 99.5% ES:", es)

# Step 7: Plot the loss distribution histogram, GPD fitted tail, and theoretical t pdf

plt.figure(figsize=(12, 8))

# Plot histogram of losses in gray
plt.hist(losses, bins=50, density=True, alpha=0.6, color='gray', label='Loss Distribution')

# Plot GPD fitted tail for x >= threshold (red line)
x = np.linspace(threshold, losses.max(), 100)
tail_pdf = genpareto.pdf(x - threshold, xi, loc=0, scale=beta)
plt.plot(x, tail_pdf * (n_exceed/n), 'r-', lw=2, label='GPD Fitted Tail')

# Plot vertical line for VaR
plt.axvline(var, color='k', linestyle='--', label='99.5% VaR')

# Plot the theoretical t-distribution pdf for losses in gray
# Given returns ~ t(df=3)/100, and losses = -returns,
# the pdf for losses: f_loss(l) = 100 * t.pdf(100 * l, df=3)
x_theory = np.linspace(losses.min(), losses.max(), 200)
pdf_t = 100 * t.pdf(100 * x_theory, df=3)
plt.plot(x_theory, pdf_t, color='gray', lw=2, label='Theoretical t PDF')

# Annotate the plot with key values; use LaTeX formatting for parameters.
textstr = '\n'.join((
    r'$u$ (95th percentile): {:.4f}'.format(threshold),
    'Number of exceedances: {}'.format(n_exceed),
    r'GPD Shape ($\xi$): {:.4f}'.format(xi),
    r'GPD Scale ($\beta$): {:.4f}'.format(beta),
    r'Estimated 99.5% VaR: {:.4f}'.format(var),
    r'Estimated 99.5% ES: {:.4f}'.format(es)
))
# 使用灰色背景 (lightgray) 來顯示文字框
props = dict(boxstyle='round', facecolor='lightgray', alpha=0.5)
plt.text(0.65, 0.75, textstr, transform=plt.gca().transAxes, fontsize=12,
         verticalalignment='top', bbox=props)

plt.xlabel("Loss")
plt.ylabel("Density")
plt.title("EVT (POT Method): GPD Fitting, Tail-Related Risk Estimation & Theoretical t PDF")
plt.legend()
plt.show()

```

![](https://cdn.hashnode.com/res/hashnode/image/upload/v1740319139388/6ea0e32e-bba6-4970-99a1-b73a238daeb2.png align="center")

圖中可見：

- **POT 方法（GPD 擬合）**：  
  - 基於超過 95% 閾值的數據進行擬合，利用有限的極端數據來外推尾部風險。  
  - 參數估計（例如形狀參數 \\(\xi\\)）可能因極端數據較少而導致尾部估計更重，結果 VaR 與 ES 值較高。  
  - 這種方法更保守，特別適合捕捉實際極端事件下可能出現的嚴重損失。

- **t 分佈理論估計**：  
  - 假定數據完全服從已知的 t 分佈（本例中 df=3, scale=1/100），直接利用分佈的分位數與密度函數計算 VaR 和 ES。  
  - 該方法依賴於假設分佈的尾部行為，通常較為樂觀，結果 VaR 與 ES 值較低。

POT 方法能捕捉到尾部數據中的極端現象，對於風險管理來說提供了更保守的估計；而 t 分佈假設可能低估了極端事件的風險。兩者各有合理性，具體選擇應根據數據特性和投資者對極端損失風險的敏感度而定。

---

# 5. 結論

1. **風險指標定義：**  
   VaR 為在給定置信水平下預計不會超過的最大損失，而 ES 為損失超過 VaR 時的平均損失。這兩個指標是極端風險管理中重要的量化工具。

2. **EVT 的作用：**  
   EVT 為捕捉極端（尾部）風險提供了理論基礎，不依賴於整體分佈假設。兩種主要方法中，Block Maxima 方法通過劃分區塊提取極值，而 POT 方法則利用所有超過高閾值的數據進行建模。

3. **POT 方法推導與外推：**  
   以 POT 方法為例，我們定義超額數據 \\(Y=X-u\\) 並假設其條件分佈滿足正則變化性質，從而推導出廣義帕累托分佈（GPD）的形式：
   
   \\[
   F_u(y) \approx 1-\left(1+\frac{\xi y}{\beta}\right)^{-1/\xi}
   \\]
   再根據該模型，利用超過閾值的數據數量 \\(n_u\\) 與總樣本數 \\(n\\)，推導出：
   - **VaR 的公式：**
     
     \\[
     VaR_\alpha = u + \frac{\beta}{\xi}\left[\left(\frac{n}{n_u}(1-\alpha)\right)^{-\xi}-1\right]
     \\]
   - **ES 的公式：**
     
     \\[
     ES_\alpha = u + \frac{\beta}{1-\xi}\left[\left(\frac{n}{n_u}(1-\alpha)\right)^{-\xi}-1\right]
     \\]

4. **參數估計：**  
   POT 方法中常用 MLE 估計 GPD 的參數 \\(\xi\\) 與 \\(\beta\\)；而 Block Maxima 方法則直接利用區塊極值進行 GEV 擬合。實務中，由於 POT 方法數據利用率高且實現簡單，常被選用。