# ===================================================================
# MODULE 5, LESSON 4: CREDIT RISK MODELLING – PD, LGD, EAD
# ===================================================================

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from scipy.stats import norm
import warnings
warnings.filterwarnings('ignore')

# Set style
sns.set_style("whitegrid")
np.random.seed(42)

print("="*70)
print("CREDIT RISK MODELLING – PD, LGD, EAD AND IFRS 9")
print("="*70)

# ----------------------------------------------------------------
# PART A: GENERATE LOAN PORTFOLIO DATA
# ----------------------------------------------------------------

n_loans = 10000

# Borrower characteristics
income = np.random.gamma(4, 20, n_loans) + 30  # $000s
dti = np.random.beta(2, 4, n_loans) * 50
credit_score = np.random.normal(700, 50, n_loans).clip(550, 850)
loan_amount = np.random.gamma(4, 40, n_loans) + 30  # $000s
term = np.random.choice([12, 24, 36, 48, 60, 72], n_loans, p=[0.1, 0.1, 0.2, 0.2, 0.2, 0.2])

# Collateral type (0 = unsecured, 1 = secured)
collateral = np.random.binomial(1, 0.6, n_loans)

# Generate PD using a logistic model
log_odds_pd = -4.0 + 0.05*dti - 0.006*credit_score + 0.005*(loan_amount/10) - 0.3*collateral
pd_true = 1/(1+np.exp(-log_odds_pd))
default = np.random.binomial(1, pd_true)

# Generate LGD (depends on collateral, recovery environment, and economic conditions)
# LGD = 1 - Recovery Rate
lgd_base = np.where(collateral == 1, 0.35, 0.65)  # secured lower LGD
lgd_noise = np.random.normal(0, 0.1, n_loans).clip(-0.15, 0.15)
lgd = (lgd_base + lgd_noise).clip(0.05, 0.95)

# EAD = loan amount + undrawn (if revolving)
ead = loan_amount  # For term loans, EAD = loan amount

# Create DataFrame
df = pd.DataFrame({
    'income': income,
    'dti': dti,
    'credit_score': credit_score,
    'loan_amount': loan_amount,
    'term': term,
    'collateral': collateral,
    'default': default,
    'lgd': lgd,
    'ead': ead,
    'pd_true': pd_true
})

print("Portfolio Summary:")
print(f"  Total Loans: {len(df):,}")
print(f"  Default Rate: {df['default'].mean():.2%}")
print(f"  Average PD: {df['pd_true'].mean():.2%}")
print(f"  Average LGD: {df['lgd'].mean():.2%}")
print(f"  Average EAD: ${df['ead'].mean():,.2f}")

# ----------------------------------------------------------------
# PART B: PD ESTIMATION USING LOGISTIC REGRESSION
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART B: PD Estimation (Logistic Regression)")
print("-"*60)

# Train-test split
X = df[['income', 'dti', 'credit_score', 'loan_amount', 'term', 'collateral']]
y = df['default']
from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=42)

# Standardise
from sklearn.preprocessing import StandardScaler
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# Logistic Regression
model_pd = LogisticRegression(random_state=42)
model_pd.fit(X_train_scaled, y_train)

# Predict PDs
pd_pred_train = model_pd.predict_proba(X_train_scaled)[:, 1]
pd_pred_test = model_pd.predict_proba(X_test_scaled)[:, 1]

# Evaluate
auc_pd = roc_auc_score(y_test, pd_pred_test)
print(f"PD Model AUC: {auc_pd:.4f}")

# Feature coefficients
coef_df = pd.DataFrame({
    'Feature': X.columns,
    'Coefficient': model_pd.coef_[0]
}).sort_values('Coefficient', ascending=False)
print("\nPD Model Coefficients:")
print(coef_df.to_string(index=False))

print("\nInterpretation:")
print("  • DTI and loan amount increase PD (positive coefficients).")
print("  • Credit score, income, and collateral decrease PD (negative coefficients).")

# ----------------------------------------------------------------
# PART C: LGD ESTIMATION
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART C: LGD Estimation")
print("-"*60)

# For demonstration, we use the true LGD (in practice, estimated from recovery data)

# LGD by collateral type
lgd_by_collateral = df.groupby('collateral')['lgd'].agg(['mean', 'std', 'count'])
print("\nLGD by Collateral Type:")
print(lgd_by_collateral)

# LGD segmentation
df['lgd_segment'] = pd.cut(df['lgd'], bins=[0, 0.2, 0.4, 0.6, 0.8, 1], 
                           labels=['0-20%', '20-40%', '40-60%', '60-80%', '80-100%'])

fig, ax = plt.subplots(figsize=(10, 6))
lgd_counts = df['lgd_segment'].value_counts().sort_index()
ax.bar(lgd_counts.index, lgd_counts.values, color='blue', alpha=0.7)
ax.set_xlabel('LGD Segment')
ax.set_ylabel('Number of Loans')
ax.set_title('LGD Distribution', fontsize=12)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('lgd_distribution.png', dpi=300)
plt.show()

# LGD Regression (predict LGD from features)
print("\nLGD Regression (OLS) – Predicting LGD from Borrower Characteristics")
from sklearn.linear_model import LinearRegression
model_lgd = LinearRegression()
model_lgd.fit(X_train_scaled, df.loc[X_train.index, 'lgd'])
lgd_pred = model_lgd.predict(X_test_scaled)

# Evaluate
from sklearn.metrics import mean_absolute_error, mean_squared_error
mae_lgd = mean_absolute_error(df.loc[X_test.index, 'lgd'], lgd_pred)
rmse_lgd = np.sqrt(mean_squared_error(df.loc[X_test.index, 'lgd'], lgd_pred))
print(f"LGD Model MAE: {mae_lgd:.4f}")
print(f"LGD Model RMSE: {rmse_lgd:.4f}")

# ----------------------------------------------------------------
# PART D: EAD ESTIMATION
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART D: EAD Estimation")
print("-"*60)

# For term loans, EAD = loan amount
# For revolving, we would model drawdowns
print("For term loans, EAD = Outstanding Balance.")
print(f"Total EAD: ${df['ead'].sum():,.2f}")

# Simulate a revolving portfolio with CCF
revolving_indicator = np.random.binomial(1, 0.3, n_loans)
df['is_revolving'] = revolving_indicator.astype(bool)
df['undrawn'] = np.where(df['is_revolving'], df['ead'] * np.random.uniform(0.2, 0.8, n_loans), 0)
df['ccf'] = np.where(df['is_revolving'], np.random.uniform(0.3, 0.8, n_loans), 0.75)
df['ead_adjusted'] = df['ead'] + df['undrawn'] * df['ccf']

print("\nEAD for Revolving Facilities:")
print(df[df['is_revolving']][['ead', 'undrawn', 'ccf', 'ead_adjusted']].head(5).round(2))

# ----------------------------------------------------------------
# PART E: EXPECTED LOSS CALCULATION
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART E: Expected Loss (EL) Calculation")
print("-"*60)

# Calculate Expected Loss for each loan
df['pd'] = pd_pred_train if len(pd_pred_train) == len(df.loc[X_train.index]) else pd_pred_test
# For this demo, use the actual PD from the model on the training set
# We'll use the full dataset predictions for demonstration
df['pd_model'] = model_pd.predict_proba(scaler.transform(X))[:, 1]
df['el'] = df['pd_model'] * df['lgd'] * df['ead_adjusted']

print(f"Total Expected Loss: ${df['el'].sum():,.2f}")
print(f"Average EL per Loan: ${df['el'].mean():,.2f}")

# Aggregate by segment
el_by_collateral = df.groupby('collateral').agg({
    'el': 'sum',
    'ead_adjusted': 'sum',
    'pd_model': 'mean',
    'lgd': 'mean'
})
el_by_collateral['el_rate'] = el_by_collateral['el'] / el_by_collateral['ead_adjusted']
print("\nEL by Collateral Type:")
print(el_by_collateral.round(4))

# ----------------------------------------------------------------
# PART F: ECONOMIC CAPITAL – UNEXPECTED LOSS (VASICEK MODEL)
# -----------------------------------------------------------------

print("\n" + "-"*60)
print("PART F: Economic Capital (Unexpected Loss)")
print("-"*60)

def vasicek_capital(pd, lgd, ead, rho=0.12, confidence=0.999):
    """
    Calculate regulatory capital using the Vasicek model.
    """
    z = norm.ppf(confidence)
    pd_adj = norm.cdf((norm.ppf(pd) + np.sqrt(rho) * z) / np.sqrt(1 - rho))
    capital = lgd * (pd_adj - pd) * ead
    return capital

# Apply to each loan
df['ul'] = vasicek_capital(df['pd_model'], df['lgd'], df['ead_adjusted'])
df['ul'] = np.maximum(df['ul'], 0)  # No negative capital

total_ul = df['ul'].sum()
print(f"Total Unexpected Loss (99.9% confidence): ${total_ul:,.2f}")
print(f"Economic Capital Ratio: {total_ul / df['ead_adjusted'].sum():.2%}")

# ----------------------------------------------------------------
# PART G: IFRS 9 – ECL CALCULATION (STAGE 1 vs STAGE 2)
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART G: IFRS 9 – ECL Calculation")
print("-"*60)

# Define SICR threshold (significant increase in credit risk)
sicr_threshold = 0.5  # PD increase of 50% or more from origination
df['pd_origination'] = df['pd_model'] * np.random.uniform(0.5, 0.9, n_loans)  # Simulated origination PD
df['pd_increase'] = (df['pd_model'] - df['pd_origination']) / df['pd_origination']
df['stage'] = 1  # Default Stage 1
df.loc[df['pd_increase'] > sicr_threshold, 'stage'] = 2  # Stage 2
df.loc[df['default'] == 1, 'stage'] = 3  # Stage 3 (defaulted)

# Lifetime PD for Stage 2/3 (approximation: multiply 12-month PD by average life)
avg_loan_life = 5  # years
df['pd_lifetime'] = df['pd_model'] * avg_loan_life
df['pd_lifetime'] = np.minimum(df['pd_lifetime'], 1)  # Cap at 1

# ECL calculation by stage
def calculate_ecl(row):
    if row['stage'] == 1:
        return row['pd_model'] * row['lgd'] * row['ead_adjusted']  # 12-month ECL
    else:  # Stages 2 and 3
        return row['pd_lifetime'] * row['lgd'] * row['ead_adjusted']  # Lifetime ECL

df['ecl_ifrs9'] = df.apply(calculate_ecl, axis=1)

print("\nIFRS 9 ECL by Stage:")
stage_summary = df.groupby('stage').agg({
    'ecl_ifrs9': 'sum',
    'ead_adjusted': 'sum',
    'pd_model': 'mean',
    'count': 'size'
})
stage_summary['ecl_rate'] = stage_summary['ecl_ifrs9'] / stage_summary['ead_adjusted']
print(stage_summary.round(4))

print(f"\nTotal IFRS 9 ECL Provision: ${df['ecl_ifrs9'].sum():,.2f}")

# Compare with CECL (lifetime ECL for all loans)
df['ecl_cecl'] = df['pd_lifetime'] * df['lgd'] * df['ead_adjusted']
print(f"Total CECL ECL Provision: ${df['ecl_cecl'].sum():,.2f}")
print(f"Difference (CECL - IFRS 9): ${(df['ecl_cecl'].sum() - df['ecl_ifrs9'].sum()):,.2f}")

# ----------------------------------------------------------------
# PART H: VISUALISATION AND BUSINESS INSIGHTS
# ----------------------------------------------------------------

fig, axes = plt.subplots(2, 2, figsize=(14, 12))

# PD Distribution
ax = axes[0, 0]
ax.hist(df['pd_model'], bins=50, alpha=0.7, edgecolor='black', color='blue')
ax.axvline(df['pd_model'].mean(), color='red', linestyle='--', label=f"Mean PD: {df['pd_model'].mean():.2%}")
ax.set_xlabel('Probability of Default (PD)')
ax.set_ylabel('Frequency')
ax.set_title('PD Distribution')
ax.legend()
ax.grid(True, alpha=0.3)

# EL by Stage
ax = axes[0, 1]
stage_el = df.groupby('stage')['ecl_ifrs9'].sum()
ax.bar(stage_el.index, stage_el.values, color=['green', 'orange', 'red'])
ax.set_xlabel('IFRS 9 Stage')
ax.set_ylabel('Total ECL ($)')
ax.set_title('ECL by Stage')
ax.grid(True, alpha=0.3)

# LGD vs PD Scatter
ax = axes[1, 0]
scatter = ax.scatter(df['pd_model'], df['lgd'], c=df['ead_adjusted'], cmap='viridis', alpha=0.5, s=10)
ax.set_xlabel('PD')
ax.set_ylabel('LGD')
ax.set_title('PD vs LGD (colored by EAD)')
plt.colorbar(scatter, ax=ax, label='EAD ($)')
ax.grid(True, alpha=0.3)

# EL vs EAD
ax = axes[1, 1]
ax.scatter(df['ead_adjusted'], df['el'], alpha=0.3, s=5, color='blue')
ax.set_xlabel('EAD ($)')
ax.set_ylabel('Expected Loss ($)')
ax.set_title('EL vs EAD')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.savefig('credit_risk_analysis.png', dpi=300)
plt.show()

# ----------------------------------------------------------------
# PART I: SUMMARY REPORT
# ----------------------------------------------------------------

print("\n" + "="*70)
print("PART I: CREDIT RISK SUMMARY REPORT")
print("="*70)

report = pd.DataFrame({
    'Metric': [
        'Number of Loans',
        'Total Exposure (EAD)',
        'Average PD',
        'Average LGD',
        'Expected Loss (EL)',
        'EL / EAD Ratio',
        'Unexpected Loss (UL)',
        'UL / EAD Ratio',
        'IFRS 9 ECL Provision',
        'CECL Provision',
        'Stage 1 Loans',
        'Stage 2 Loans',
        'Stage 3 Loans'
    ],
    'Value': [
        f"{len(df):,}",
        f"${df['ead_adjusted'].sum():,.2f}",
        f"{df['pd_model'].mean():.2%}",
        f"{df['lgd'].mean():.2%}",
        f"${df['el'].sum():,.2f}",
        f"{df['el'].sum() / df['ead_adjusted'].sum():.2%}",
        f"${df['ul'].sum():,.2f}",
        f"{df['ul'].sum() / df['ead_adjusted'].sum():.2%}",
        f"${df['ecl_ifrs9'].sum():,.2f}",
        f"${df['ecl_cecl'].sum():,.2f}",
        f"{(df['stage'] == 1).sum():,}",
        f"{(df['stage'] == 2).sum():,}",
        f"{(df['stage'] == 3).sum():,}"
    ]
})
print(report.to_string(index=False))