Python 使用 FFT 逼近聚合损失随机变量的 CDF

下面是我几周前分配给我的一个类的python代码，我一直无法成功调试.问题在于如何使用FFT找到总损失随机变量的风险值(即p%分位数).我们给出了一个清晰的数学过程，通过这个过程我们可以得到总损失随机变量的离散化CDF的估计.然而，我的结果非常糟糕，我犯了一些错误，即使在调试了几个小时的代码之后，我也找不到这些错误.

给出了总损失随机变量S，其中S=sum(X_i for i in range(N))与r=5, beta=.2呈负二项分布，X_i与theta=1呈指数分布.这个参数化的概率母函数是P(z)=[1-\beta(z-1)]^{-r}.

我们被要求估计S人的分布情况

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

得到的向量应该是S的每个这样的区间的概率质量的近似值.我从以前的方法中知道，95%的VaR应该是~4，99.9%的VaR应该是~10.但我的代码返回的结果毫无意义.一般来说，我的ECDF达到的指数&gt；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