q = 0.999
p = 0.001
e = 1.0
number_samples = 10000000
def capital_analytic(p, e, q=0.999):
return e * (np.heaviside(p + q -1 , 1) - p)
def capital_MC(p, e, number_samples, q=0.999):
rng = np.random.default_rng()
loss_indicator = np.heaviside(p-rng.random(number_samples), 1)
loss = e * loss_indicator
quantile = np.quantile(loss, q)
expected_loss = np.mean(loss)
return quantile - expected_loss
capital_simulation = capital_MC(p, e, number_samples, q)
capital_exact = capital_analytic(p, e, q)
print(f"Capital: MC={capital_simulation}, analytic={capital_exact}, rel. diff.: {(capital_simulation - capital_exact)/capital_exact:.2%}")
ps = np.linspace(0.0, 0.01, 1001)
c_sim = [capital_MC(p, e, 10000, q) for p in ps]
c_exact = [capital_analytic(p, e, q) for p in ps]
#plt.plot(ps, c_sim, c="red")
plt.plot(ps, c_exact, c="blue", alpha=0.4)
plt.xlabel("p")
plt.ylabel("capital/EAD")
plt.title("Capital for Portfolio with Single Position")
#plt.vlines(1-q, -0.1, 1.1, colors="r", linestyles="dotted", label=f"p={1-q}")
plt.axvline(1-q, color="red", linestyle=":")
plt.text(1-q,0.5, f"p=1-q", color="red")
plt.show()
