# ===================================================================
# MODULE 4, LESSON 7: CREDIT RISK MODELLING WITH AI
# ===================================================================

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from xgboost import XGBClassifier
from sklearn.metrics import (roc_auc_score, classification_report, confusion_matrix,
                             roc_curve, mean_squared_error, r2_score)
from scipy.stats import norm
import warnings
warnings.filterwarnings('ignore')

print("="*70)
print("CREDIT RISK MODELLING WITH AI")
print("="*70)

# ----------------------------------------------------------------
# PART A: GENERATE CREDIT DATA
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART A: Generating Credit Data")
print("-"*60)

np.random.seed(42)
n_loans = 10000

# Generate borrower features
credit_data = pd.DataFrame({
    'loan_id': range(1, n_loans + 1),
    'age': np.random.normal(45, 15, n_loans).clip(18, 80).astype(int),
    'income': np.random.gamma(5, 20, n_loans) + 20,
    'credit_score': np.random.normal(700, 50, n_loans).clip(550, 850).astype(int),
    'dti': np.random.beta(2, 5, n_loans) * 60,
    'loan_amount': np.random.gamma(4, 50, n_loans) + 30,
    'loan_term': np.random.choice([12, 24, 36, 48, 60, 72], n_loans),
    'employment_years': np.random.gamma(3, 5, n_loans).clip(0, 30).astype(int),
    'home_owner': np.random.binomial(1, 0.65, n_loans),
    'marital_status': np.random.choice([0, 1, 2], n_loans, p=[0.35, 0.45, 0.20]),
    'education': np.random.choice([0, 1, 2, 3], n_loans, p=[0.15, 0.25, 0.35, 0.25]),
    'previous_defaults': np.random.choice([0, 1, 2], n_loans, p=[0.8, 0.15, 0.05])
})

# Generate default based on features
log_odds = (-4.5 + 0.04 * credit_data['dti'] - 0.005 * credit_data['credit_score']
            + 0.01 * (credit_data['loan_amount']/1000) + 0.02 * credit_data['employment_years']
            - 0.01 * credit_data['age'] + 0.3 * credit_data['home_owner']
            - 0.5 * credit_data['marital_status'] - 0.3 * credit_data['education']
            + 0.2 * credit_data['previous_defaults'])

# Add non-linearity
log_odds += 0.0003 * credit_data['dti']**2
log_odds -= 0.00001 * credit_data['credit_score'] * credit_data['dti']

prob_default = 1 / (1 + np.exp(-log_odds))
credit_data['default'] = np.random.binomial(1, prob_default)

# Generate LGD based on home ownership
credit_data['lgd'] = np.where(credit_data['home_owner'] == 1,
                             np.random.normal(0.35, 0.08, n_loans).clip(0.10, 0.70),
                             np.random.normal(0.65, 0.10, n_loans).clip(0.20, 0.90))

# EAD = loan amount
credit_data['ead'] = credit_data['loan_amount']

print(f"Generated {len(credit_data)} loan applications")
print(f"Default rate: {credit_data['default'].mean():.2%}")
print(f"Average LGD: {credit_data['lgd'].mean():.2%}")

# ----------------------------------------------------------------
# PART B: EXPLORATORY ANALYSIS
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART B: Exploratory Analysis")
print("-"*60)

# Default rate by key features
default_summary = credit_data.groupby('default').agg({
    'credit_score': 'mean',
    'dti': 'mean',
    'income': 'mean',
    'loan_amount': 'mean'
}).round(2)

print("Default Summary:")
print(default_summary)

# Visualise
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# Credit Score vs Default
ax = axes[0, 0]
sns.boxplot(data=credit_data, x='default', y='credit_score', ax=ax)
ax.set_xticklabels(['No Default', 'Default'])
ax.set_title('Credit Score by Default')

# DTI vs Default
ax = axes[0, 1]
sns.boxplot(data=credit_data, x='default', y='dti', ax=ax)
ax.set_xticklabels(['No Default', 'Default'])
ax.set_title('DTI by Default')

# Income vs Loan Amount
ax = axes[1, 0]
scatter = ax.scatter(credit_data['income'], credit_data['loan_amount'],
                     c=credit_data['default'], cmap='RdYlGn', alpha=0.5, s=10)
ax.set_xlabel('Income ($000s)')
ax.set_ylabel('Loan Amount ($000s)')
ax.set_title('Income vs Loan Amount')
plt.colorbar(scatter, ax=ax, label='Default')

# Correlation Heatmap
ax = axes[1, 1]
numeric_cols = ['age', 'income', 'credit_score', 'dti', 'loan_amount',
                'loan_term', 'employment_years', 'default']
corr = credit_data[numeric_cols].corr()
sns.heatmap(corr, annot=True, fmt='.2f', cmap='coolwarm', ax=ax)
ax.set_title('Correlation Matrix')

plt.tight_layout()
plt.savefig('credit_risk_eda.png', dpi=300, bbox_inches='tight')
plt.show()
print("Credit risk EDA visualisation saved as 'credit_risk_eda.png'")

# ----------------------------------------------------------------
# PART C: FEATURE ENGINEERING
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART C: Feature Engineering")
print("-"*60)

# Create engineered features
credit_data['loan_to_income'] = credit_data['loan_amount'] / credit_data['income']
credit_data['dti_credit_interaction'] = credit_data['dti'] * credit_data['credit_score'] / 1000
credit_data['dti_squared'] = credit_data['dti'] ** 2
credit_data['credit_score_category'] = pd.cut(credit_data['credit_score'],
                                             bins=[550, 600, 650, 700, 750, 850],
                                             labels=['Poor', 'Fair', 'Good', 'Very Good', 'Excellent'])

# Features for modelling
features = ['age', 'income', 'credit_score', 'dti', 'loan_amount', 'loan_term',
            'employment_years', 'home_owner', 'marital_status', 'education',
            'previous_defaults', 'loan_to_income', 'dti_credit_interaction',
            'dti_squared']

X = credit_data[features]
y = credit_data['default']

# Encode categorical features
X_encoded = pd.get_dummies(X, columns=['marital_status', 'education'], drop_first=True)

print(f"Features: {len(X_encoded.columns)}")

# ----------------------------------------------------------------
# PART D: MODEL TRAINING AND EVALUATION
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART D: Model Training and Evaluation")
print("-"*60)

# Train-test split
X_train, X_test, y_train, y_test = train_test_split(X_encoded, y, test_size=0.3, random_state=42)

# Scale features
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

# 1. Logistic Regression (Baseline)
lr_model = LogisticRegression(class_weight='balanced', max_iter=1000, random_state=42)
lr_model.fit(X_train_scaled, y_train)
y_pred_lr_proba = lr_model.predict_proba(X_test_scaled)[:, 1]
auc_lr = roc_auc_score(y_test, y_pred_lr_proba)

# 2. Random Forest
rf_model = RandomForestClassifier(n_estimators=100, max_depth=10, class_weight='balanced', random_state=42)
rf_model.fit(X_train_scaled, y_train)
y_pred_rf_proba = rf_model.predict_proba(X_test_scaled)[:, 1]
auc_rf = roc_auc_score(y_test, y_pred_rf_proba)

# 3. XGBoost
scale_pos_weight = len(y_train[y_train==0]) / len(y_train[y_train==1])
xgb_model = XGBClassifier(n_estimators=100, max_depth=6, learning_rate=0.1,
                         scale_pos_weight=scale_pos_weight,
                         random_state=42, use_label_encoder=False, eval_metric='logloss')
xgb_model.fit(X_train_scaled, y_train)
y_pred_xgb_proba = xgb_model.predict_proba(X_test_scaled)[:, 1]
auc_xgb = roc_auc_score(y_test, y_pred_xgb_proba)

print("Model Performance (AUC):")
print(f"  Logistic Regression: {auc_lr:.4f}")
print(f"  Random Forest: {auc_rf:.4f}")
print(f"  XGBoost: {auc_xgb:.4f}")

# Best model = XGBoost
y_pred = (y_pred_xgb_proba >= 0.5).astype(int)
print("\nXGBoost Classification Report:")
print(classification_report(y_test, y_pred, target_names=['No Default', 'Default']))

# Confusion Matrix
cm = confusion_matrix(y_test, y_pred)
print("\nConfusion Matrix:")
print(pd.DataFrame(cm, columns=['Pred No Default', 'Pred Default'],
                   index=['Actual No Default', 'Actual Default']))

# Feature Importance
importance_xgb = pd.DataFrame({
    'Feature': X_encoded.columns,
    'Importance': xgb_model.feature_importances_
}).sort_values('Importance', ascending=False)

print("\nTop 10 Credit Risk Predictors:")
print(importance_xgb.head(10).to_string(index=False))

# ----------------------------------------------------------------
# PART E: CALIBRATION AND STRESS TESTING
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART E: Calibration and Stress Testing")
print("-"*60)

# Calibration plot
def calibration_plot(y_true, y_pred, n_bins=10):
    """Create a calibration plot."""
    df = pd.DataFrame({'y_true': y_true, 'y_pred': y_pred})
    df['bin'] = pd.qcut(df['y_pred'], q=n_bins, labels=False, duplicates='drop')
    cal_data = df.groupby('bin').agg(
        observed=('y_true', 'mean'),
        predicted=('y_pred', 'mean'),
        count=('y_true', 'count')
    ).reset_index()
    return cal_data

cal_data = calibration_plot(y_test, y_pred_xgb_proba)

fig, ax = plt.subplots(figsize=(8, 6))
ax.plot(cal_data['predicted'], cal_data['observed'], 'bo-', linewidth=2, markersize=8, label='Model')
ax.plot([0, 1], [0, 1], 'r--', linewidth=1, label='Perfect Calibration')
ax.set_xlabel('Mean Predicted Probability')
ax.set_ylabel('Observed Default Rate')
ax.set_title('Calibration Plot - XGBoost')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('credit_calibration.png', dpi=300, bbox_inches='tight')
plt.show()
print("Calibration plot saved as 'credit_calibration.png'")

# Expected Loss calculation
credit_data['pd'] = xgb_model.predict_proba(scaler.transform(X_encoded))[:, 1]
credit_data['el'] = credit_data['pd'] * credit_data['lgd'] * credit_data['ead']

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

# Stress testing (simulated)
def stress_test(df, gdp_shock, unemp_shock):
    """Simulate stress testing on credit portfolio."""
    df_stressed = df.copy()
    # Apply macroeconomic shocks
    log_odds = np.log(df_stressed['pd'] / (1 - df_stressed['pd']))
    log_odds += -1.5 * gdp_shock + 3.0 * unemp_shock
    df_stressed['pd_stressed'] = 1 / (1 + np.exp(-log_odds))
    df_stressed['pd_stressed'] = df_stressed['pd_stressed'].clip(0.001, 0.999)
    df_stressed['el_stressed'] = df_stressed['pd_stressed'] * df_stressed['lgd'] * df_stressed['ead']
    return df_stressed

scenarios = {
    'Baseline': {'gdp': 0.0, 'unemp': 0.0},
    'Adverse': {'gdp': -0.02, 'unemp': 0.03},
    'Severe': {'gdp': -0.05, 'unemp': 0.06}
}

stress_results = []
for name, shocks in scenarios.items():
    df_stressed = stress_test(credit_data, shocks['gdp'], shocks['unemp'])
    el_stressed = df_stressed['el_stressed'].sum()
    stress_results.append({
        'Scenario': name,
        'EL': el_stressed,
        'Increase %': (el_stressed / credit_data['el'].sum() - 1) * 100
    })

stress_df = pd.DataFrame(stress_results)
print("\nStress Test Results:")
print(stress_df.to_string(index=False))

# ----------------------------------------------------------------
# PART F: REGULATORY METRICS
# ----------------------------------------------------------------

print("\n" + "-"*60)
print("PART F: Regulatory Metrics")
print("-"*60)

# KS Statistic
from scipy.stats import ks_2samp
scores_good = y_pred_xgb_proba[y_test == 0]
scores_bad = y_pred_xgb_proba[y_test == 1]
ks_stat, _ = ks_2samp(scores_good, scores_bad)

# Gini Coefficient
gini = 2 * auc_xgb - 1

print(f"KS Statistic: {ks_stat:.4f}")
print(f"Gini Coefficient: {gini:.4f}")
print(f"AUC: {auc_xgb:.4f}")

# Regulatory benchmarks
regulatory_benchmarks = {
    'Metric': ['AUC', 'KS', 'Gini', 'Calibration (H-L)'],
    'Acceptable': ['> 0.75', '> 0.30', '> 0.50', 'p > 0.05'],
    'Good': ['> 0.80', '> 0.35', '> 0.60', 'p > 0.10'],
    'Excellent': ['> 0.85', '> 0.40', '> 0.70', 'p > 0.20']
}

benchmark_df = pd.DataFrame(regulatory_benchmarks)
print("\nRegulatory Benchmarks:")
print(benchmark_df.to_string(index=False))

# ----------------------------------------------------------------
# PART G: SUMMARY AND RECOMMENDATIONS
# ----------------------------------------------------------------

print("\n" + "="*70)
print("PART G: Summary and Recommendations")
print("="*70)

print("""
Credit Risk Modelling with AI – Key Takeaways:

1. Credit risk includes PD, LGD, EAD, Expected Loss, and Unexpected Loss.
2. AI models (XGBoost, Random Forest) outperform traditional logistic regression.
3. Alternative data improves credit risk assessment for underserved populations.
4. Model validation: conceptual soundness, data quality, performance, stability.
5. Regulatory framework: Basel III, IFRS 9/CECL, ECOA/Fair Lending.
6. Key metrics: AUC, KS, Gini, calibration, stress testing.
7. Explainability (SHAP) is essential for regulatory compliance.

Recommendations:
  - Use ensemble models (XGBoost) for credit scoring.
  - Incorporate alternative data for financial inclusion.
  - Regularly validate and monitor credit risk models.
  - Ensure model explainability for regulatory compliance.
  - Conduct stress testing under multiple scenarios.
  - Maintain clear documentation for regulatory submissions.
""")

print("="*70)
print("END OF LESSON 7 – MODULE 4")
print("="*70)