import numpy as np
from scipy.stats import norm
# -----------------------------
# PARAMETERS
# -----------------------------
N = 125 # number of names
rho = 0.2 # asset correlation
p_default = 0.02 # annual default probability
recovery = 0.4
tranches = [(0.0, 0.03), (0.03, 0.07), (0.07, 1.0)]
n_sims = 20000
# -----------------------------
# PRECOMPUTE DEFAULT THRESHOLD
# -----------------------------
# Step 3 concept: convert default probability into latent threshold
default_threshold = norm.ppf(p_default)
# equal weights portfolio
weights = np.ones(N) / N
# -----------------------------
# MONTE CARLO SIMULATION
# -----------------------------
tranche_losses = np.zeros((n_sims, len(tranches)))
for s in range(n_sims):
# ============================================================
# STEP 1: GENERATE SYSTEMATIC MARKET SHOCK
# This represents the shared macroeconomic factor (M)
# that drives correlated defaults across all names.
# ============================================================
M = np.random.normal()
# ============================================================
# STEP 2: GENERATE IDIOSYNCRATIC SHOCKS
# Each obligor has its own independent noise term (epsilon)
# capturing firm-specific risk not explained by the market.
# ============================================================
eps = np.random.normal(size=N)
# ============================================================
# STEP 3: BUILD CORRELATED LATENT VARIABLES (GAUSSIAN COPULA)
# Z_i = sqrt(rho)*M + sqrt(1-rho)*epsilon_i
# This introduces dependence between obligors via the
# shared market factor while preserving normal marginals.
# ============================================================
Z = np.sqrt(rho)*M + np.sqrt(1-rho)*eps
# ============================================================
# STEP 4: MAP LATENT VARIABLE TO DEFAULT EVENT
# If Z_i < Phi^{-1}(p), the obligor defaults.
# This converts continuous risk factors into binary defaults.
# ============================================================
defaults = Z < default_threshold
# -----------------------------
# PORTFOLIO LOSS CALCULATION
# -----------------------------
losses = defaults * weights * (1 - recovery)
total_loss = np.sum(losses)
# -----------------------------
# TRANCHE ALLOCATION
# -----------------------------
for t, (A, D) in enumerate(tranches):
tranche_loss = min(max(total_loss - A, 0), D - A)
tranche_losses[s, t] = tranche_loss
# -----------------------------
# RESULTS
# -----------------------------
expected_tranche_loss = tranche_losses.mean(axis=0)
for i, (A, D) in enumerate(tranches):
print(f"Tranche {A*100:.1f}-{D*100:.1f}% Expected Loss: {expected_tranche_loss[i]:.4f}")
