python - 使用 FFT 来近似聚合损失随机变量的 CDF

Question

您将在下面找到我几周前收到的类作业的 Python 代码，但我无法成功调试。问题是关于使用 FFT 找到聚合损失随机变量的风险值(即 p% 分位数)。我们得到了一个清晰的数学程序，通过它我们可以获得聚合损失随机变量的离散化 CDF 的估计。然而，我的结果严重偏离，我犯了某种错误，即使在调试我的代码数小时后我也找不到。

总损失随机变量 S 是 S=sum(X_i for i in range(N))，其中 N服从 r=5, beta=.2 的负二项分布，X_i 服从 theta=1 的指数分布。此参数化的概率生成函数为 P(z)=[1-\beta(z-1)]^{-r}。

我们被要求通过

来近似 S的分布

1. 选择网格宽度 h 和整数 n 使得 r=2^n 是要离散化的元素数量 X 上，
2. 离散化X并计算处于宽度h等间隔区间的概率，
3. 将 FFT 应用于离散化的 X，
4. 将 N 的 PGF 应用于傅里叶变换的 X 的元素，
5. 对该向量应用逆 FFT。

生成的向量应该是 S 的每个此类区间的概率质量的近似值。我从以前的方法中知道 95% VaR 应该是 ~4 而 99.9% VaR 应该是~10。但是我的代码返回了无意义的结果。一般来说，我的 ECDF 达到水平 > 0.95 的指数已经太晚了，即使经过数小时的调试，我也没有找到我哪里出错了。

我也在 math stackexchange 上问过这个问题，因为这个问题在很大程度上是编程和数学的交叉点，我现在不知道这个问题是在实现方面还是我正在申请数学思想错误。

import numpy as np
from scipy.stats import expon
from scipy.fft import fft, ifft

r, beta, theta = 5, .2, 1
var_levels = [.95, .999]


def discretize_X(h: float, m: int):
    X = expon(scale=theta)
    f_X = [X.cdf(h / 2),
           *[X.cdf(j * h + h / 2) - X.cdf(j * h - h / 2) for j in range(1, m - 1)],
           X.sf((m - 1) * h - h / 2)]
    return f_X


# Probability generating function of N ~ NB(r, beta)
def PGF(z: [float, complex]):
    return (1 - beta * (z - 1)) ** (-r)


h = 1e-2
n = 10
r = 2 ** n

VaRs, TVaRs = [], []

# discretize X with (r-1) cells of width h and one final cell with the survival function at h*(r-1)
f_X = discretize_X(h, r)
phi_vec = fft(f_X)
f_tilde_vec_fft = np.array([PGF(phi) for phi in phi_vec])
f_S = np.real(ifft(f_tilde_vec_fft))
ecdf_S = np.cumsum(f_S)  # calc cumsum to get ECDF

for p in var_levels:
    var_idx = np.where(ecdf_S >= p)[0][0]  # get lowest index where ecdf_S >= p
    print("p =", p, "\nVaR idx:", var_idx)
    var = h * var_idx  # VaR should be this index times the cell width
    print("VaR:", var)
    tvar = 1 / (1 - p) * np.sum(f_S[var_idx:] * np.array([i * h for i in range(var_idx, r)]))  # TVaR should be each cell's probability times the value inside that cell

    VaRs.append(var)
    TVaRs.append(tvar)

return VaRs, TVaRs

Answer 1

不确定数学，但在代码片段中变量 r 被覆盖，并且在计算 f_tilde_vec_fft 函数时 PGF 使用 not 5 符合 r 的预期，但 1024。修复 -- 在超参数定义中将名称 r 更改为 r_nb:

r_nb, beta, theta = 5, .2, 1

还有函数 PGF:

返回 (1 - beta * (z - 1)) ** (-r_nb)

运行后，对于 VaRs，其他参数保持不变(例如 h、n 等)，我得到 [4.05, 9.06]