Consider adding credit risk model for a single portfolio example

[dafeda]

Based on section 20.2.1 of Risk Analysis book by Vose (2008).

Note that ModelRisk is able to do this kind of thing fast using Fourier transforms.

import numpy as np
from probabilit.modeling import Distribution
from probabilit.distributions import Lognormal

rng = np.random.default_rng(42)

# Portfolio parameters
n_obligors = 2135
default_prob = 0.083  # 8.30%

# Create the aggregate loss model
total_loss = 0
for i in range(n_obligors):
    # Each obligor either defaults or doesn't
    default_indicator = Distribution("bernoulli", p=default_prob)
    individual_loss = Lognormal(55, 12)
    
    # Obligor's contribution = default_indicator * individual_loss
    obligor_loss = default_indicator * individual_loss
    total_loss += obligor_loss

print(f"Number of nodes in graph: {total_loss.num_distribution_nodes()}")

print("Sampling... this will take some time due to the large computational graph")
samples = total_loss.sample(100000, random_state=rng)

def calculate_var(samples, confidence_level=0.95):
    "95% of the time, losses won't exceed this amount"
    var = np.percentile(samples, confidence_level * 100)
    return var

var_95 = calculate_var(samples)

print(f"\nResults:")
print(f"Mean Loss: {np.mean(samples):.2f}")
print(f"VaR (95%): {var_95:.2f}")
print(f"Standard deviation: {np.std(samples):.2f}")
print(f"\nExpected theoretical mean: {n_obligors * default_prob * 55:.2f}")