Approximating CDF for an aggregate loss model using FFT

Question

For a class assignment, I was asked to calculate the VaR at levels 95% and 99.9% for an aggregate loss random variable $S=\sum_{i=1}^{N}{X_i}$, where $N \sim \mathrm{NB}(r=5,\beta=\frac{1}{5})$, and $X_i \overset{iid}{\sim} \exp{(\theta=1)}$. We were to do this by using an FFT on a discretization of the continuous RV $X$.

The way to do this is by

1. choosing a grid width $h$ and an integer $n$ such that $r:=2^n$ is the number of elements to discretize $X$ on,
2. discretizing $X$ and calculating the probabilities of being in equally spaced intervals of width $h$,
3. applying the FFT to the discretized $X$,
4. applying the PGF of $N$ to the elements of the Fourier-transformed $X$,
5. applying the inverse FFT to this vector.

The resulting vector should be an approximation for the probability masses of each such interval for $S$. Below I attach my Python code. I know from previous methods that the 95% VaR ought to be ~4 and the 99.9% VaR ought to be ~10. But my code returns nonsensical results. Generally speaking, my index where the ECDF reaches levels $>0.95$ is way too late, and even after hours of debugging I have not managed to find where I am going wrong.

The probability generating function for this parametrization of a negative binomial distribution is $P(z)=[1-\beta(z-1)]^{-r}$.

I hope that this is the correct place to ask this question, as I find it is mostly math- and only minorly programming related.

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
```