使用 FFT 来近似聚合损失随机变量的 CDF
Using FFT to approximate the CDF for an aggregate loss random variable
您将在下面找到我的 python 代码,用于几周前我收到的 class 作业,但我无法成功调试。问题是关于使用 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
的分布
- 选择网格宽度
h
和一个整数 n
使得 r=2^n
是要在 X
上离散化的元素数,
- 离散化
X
并计算处于宽度 h
、 等间隔区间的概率
- 将 FFT 应用于离散化
X
,
- 将
N
的 PGF 应用于傅里叶变换的元素 X
,
- 正在对该向量应用逆向 FFT。
生成的向量应该是 S
的每个此类区间的概率质量的近似值。我从以前的方法中知道 95% VaR 应该是 ~4 而 99.9% VaR 应该是~10。但是我的代码 returns 无意义的结果。一般来说,我的 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
不确定数学,但在代码段中变量 r
被覆盖,并且在计算 f_tilde_vec_fft
函数时 PGF
使用的不是 5
作为 r
,但 1024
。修复——将超参数定义中的名称 r
更改为 r_nb
:
r_nb, beta, theta = 5, .2, 1
还有函数 PGF
:
return (1 - beta * (z - 1)) ** (-r_nb)
在 运行 之后,其他参数保持不变(例如 h
、n
等)对于 VaRs
我得到 [4.05, 9.06]
您将在下面找到我的 python 代码,用于几周前我收到的 class 作业,但我无法成功调试。问题是关于使用 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
的分布
- 选择网格宽度
h
和一个整数n
使得r=2^n
是要在X
上离散化的元素数, - 离散化
X
并计算处于宽度h
、 等间隔区间的概率
- 将 FFT 应用于离散化
X
, - 将
N
的 PGF 应用于傅里叶变换的元素X
, - 正在对该向量应用逆向 FFT。
生成的向量应该是 S
的每个此类区间的概率质量的近似值。我从以前的方法中知道 95% VaR 应该是 ~4 而 99.9% VaR 应该是~10。但是我的代码 returns 无意义的结果。一般来说,我的 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
不确定数学,但在代码段中变量 r
被覆盖,并且在计算 f_tilde_vec_fft
函数时 PGF
使用的不是 5
作为 r
,但 1024
。修复——将超参数定义中的名称 r
更改为 r_nb
:
r_nb, beta, theta = 5, .2, 1
还有函数 PGF
:
return (1 - beta * (z - 1)) ** (-r_nb)
在 运行 之后,其他参数保持不变(例如 h
、n
等)对于 VaRs
我得到 [4.05, 9.06]