# =============================================================================
# MÓDULO 1: ANÁLISIS EXPLORATORIO — CREDIT RISK
# Autor: Victor Miranda | Especialista en Riesgos Financieros
# Entorno: Google Colab
# =============================================================================
# !pip install imbalanced-learn  # Ejecutar en Colab si necesario

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from scipy import stats
from sklearn.experimental import enable_iterative_imputer
from sklearn.impute import IterativeImputer
from imblearn.over_sampling import SMOTE

np.random.seed(42)
sns.set_style("whitegrid")

# --- GENERACIÓN DE DATOS SINTÉTICOS (10,000 clientes) ---
n = 10000
default_rate = 0.035
n_def = int(n * default_rate)
target = np.concatenate([np.ones(n_def), np.zeros(n - n_def)])
np.random.shuffle(target)

data = pd.DataFrame()
data['default'] = target.astype(int)

# Edad: defaults más jóvenes
data['edad'] = np.where(data['default']==1,
    np.random.normal(32,8,n), np.random.normal(42,12,n)
).clip(18,75).astype(int)

# Ingreso anual (€)
data['ingreso_anual'] = np.where(data['default']==1,
    np.random.lognormal(10.0,0.5,n), np.random.lognormal(10.5,0.45,n)
).round(2)

# DTI (Debt-to-Income)
data['dti'] = np.where(data['default']==1,
    np.random.beta(5,5,n)*0.6+0.2, np.random.beta(2,5,n)*0.5+0.05
).clip(0.01,0.95).round(4)

# Antigüedad laboral
data['antiguedad_laboral'] = np.where(data['default']==1,
    np.random.exponential(3,n), np.random.exponential(7,n)
).clip(0,40).round(1)

# Número de productos
data['num_productos'] = np.where(data['default']==1,
    np.random.poisson(1.5,n), np.random.poisson(3.0,n)).clip(0,10)

# Score bureau (300-850)
data['score_bureau'] = np.where(data['default']==1,
    np.random.normal(520,80,n), np.random.normal(700,70,n)
).clip(300,850).astype(int)

# Meses desde último atraso
data['meses_ultimo_atraso'] = np.where(data['default']==1,
    np.random.exponential(6,n), np.random.exponential(24,n)
).clip(0,120).astype(int)

# Tipo empleo (categórica)
tipos = ['asalariado','autonomo','funcionario','temporal']
data['tipo_empleo'] = np.where(data['default']==1,
    np.random.choice(tipos, n, p=[.35,.25,.10,.30]),
    np.random.choice(tipos, n, p=[.45,.15,.25,.15]))

# Introducir missings y outliers
for col, pct in [('ingreso_anual',.05),('antiguedad_laboral',.08),('score_bureau',.03)]:
    data.loc[np.random.random(n) < pct, col] = np.nan
data.loc[np.random.choice(n,15,replace=False), 'ingreso_anual'] = \
    np.random.uniform(500000,2000000,15)

print(f"Shape: {data.shape} | Default rate: {data['default'].mean():.2%}")
print(f"Missing:\n{data.isnull().sum()}")

# --- ESTADÍSTICAS DESCRIPTIVAS POR GRUPO ---
vars_num = ['edad','ingreso_anual','dti','antiguedad_laboral',
            'num_productos','score_bureau','meses_ultimo_atraso']
print("\n--- NO DEFAULT ---")
print(data[data['default']==0][vars_num].describe().round(2))
print("\n--- DEFAULT ---")
print(data[data['default']==1][vars_num].describe().round(2))

# --- ANÁLISIS DE MOMENTOS ---
moments = pd.DataFrame({
    'Variable': vars_num,
    'Skewness': [data[v].skew() for v in vars_num],
    'Exc_Kurtosis': [data[v].kurtosis() for v in vars_num]
})
print("\nMOMENTOS:"); print(moments.to_string(index=False))

# --- VISUALIZACIONES ---
fig, axes = plt.subplots(2,3, figsize=(16,10))
fig.suptitle('Distribuciones Default vs No-Default', fontsize=14, fontweight='bold')
for idx, var in enumerate(['edad','ingreso_anual','dti',
                           'antiguedad_laboral','score_bureau','meses_ultimo_atraso']):
    ax = axes[idx//3, idx%3]
    for lbl, clr in [(0,'#2196F3'),(1,'#F44336')]:
        ax.hist(data[data['default']==lbl][var].dropna(), bins=40,
                alpha=0.6, color=clr, density=True,
                label='Default' if lbl else 'No Default')
    ax.set_title(var, fontweight='bold'); ax.legend(fontsize=9)
plt.tight_layout(); plt.show()

# --- IMPUTACIÓN MICE ---
imputer = IterativeImputer(max_iter=10, random_state=42)
data[vars_num] = imputer.fit_transform(data[vars_num])
print(f"\nMissing post-MICE: {data[vars_num].isnull().sum().sum()}")

# --- OUTLIERS: MAHALANOBIS ---
X = data[vars_num].values
mu = X.mean(axis=0)
cov_inv = np.linalg.inv(np.cov(X, rowvar=False))
mahal = np.array([np.sqrt((x-mu) @ cov_inv @ (x-mu).T) for x in X])
threshold = np.sqrt(stats.chi2.ppf(0.975, len(vars_num)))
print(f"Outliers Mahalanobis: {(mahal>threshold).sum()} ({(mahal>threshold).mean():.2%})")

# Winsorización percentiles 1-99
for v in vars_num:
    q1, q99 = data[v].quantile([0.01,0.99])
    data[v] = data[v].clip(q1, q99)

# --- SMOTE ---
smote = SMOTE(sampling_strategy=0.3, random_state=42, k_neighbors=5)
X_res, y_res = smote.fit_resample(data[vars_num].values, data['default'].values)
print(f"\nSMOTE: No-Def={(y_res==0).sum()}, Def={(y_res==1).sum()}")
print(f"Ratio: {(y_res==1).sum()/(y_res==0).sum():.1%}")

# =============================================================================
# MÓDULO 2: FEATURE ENGINEERING — WOE & IV
# Autor: Victor Miranda
# =============================================================================
# (Continuamos con 'data' del Módulo 1, post-imputación)

from sklearn.preprocessing import StandardScaler

# --- ESTANDARIZACIÓN ---
vars_std = ['edad','ingreso_anual','dti','antiguedad_laboral',
            'score_bureau','meses_ultimo_atraso']

scaler = StandardScaler()
data_z = pd.DataFrame(scaler.fit_transform(data[vars_std]),
                      columns=[f'{v}_z' for v in vars_std])
print("Z-Score (mean≈0, std≈1):")
print(data_z.describe().loc[['mean','std']].round(4))

# --- FUNCIÓN WOE & IV ---
def calculate_woe_iv(df, feature, target, n_bins=10):
    """Calcula WOE e IV para variable continua."""
    temp = df[[feature, target]].dropna().copy()
    temp['bin'] = pd.qcut(temp[feature], q=n_bins, duplicates='drop')
    
    g = temp.groupby('bin', observed=True).agg(
        total=(feature,'count'), n_def=(target,'sum')).reset_index()
    g['n_nd'] = g['total'] - g['n_def']
    
    tot_d, tot_nd = g['n_def'].sum(), g['n_nd'].sum()
    eps = 0.0001
    g['pct_d'] = (g['n_def']/tot_d).clip(lower=eps)
    g['pct_nd'] = (g['n_nd']/tot_nd).clip(lower=eps)
    g['woe'] = np.log(g['pct_nd'] / g['pct_d'])
    g['iv_bin'] = (g['pct_nd'] - g['pct_d']) * g['woe']
    g['def_rate'] = g['n_def'] / g['total']
    return g, g['iv_bin'].sum()

# Calcular IV para todas las variables
iv_results, woe_tables = {}, {}
for var in vars_std:
    tbl, iv = calculate_woe_iv(data, var, 'default')
    iv_results[var] = iv
    woe_tables[var] = tbl

iv_rank = pd.DataFrame(sorted(iv_results.items(), key=lambda x: x[1],
                               reverse=True), columns=['Variable','IV'])
iv_rank['Power'] = iv_rank['IV'].apply(
    lambda x: 'No pred.' if x<.02 else 'Débil' if x<.10
    else 'Medio' if x<.30 else 'Fuerte' if x<.50 else 'Sospechoso')
print("\nRANKING POR IV:")
print(iv_rank.to_string(index=False))

# --- VISUALIZACIÓN WOE ---
fig, axes = plt.subplots(2,3, figsize=(16,10))
fig.suptitle('Weight of Evidence por Variable', fontsize=14, fontweight='bold')
for idx, var in enumerate(vars_std):
    ax = axes[idx//3, idx%3]
    tbl = woe_tables[var]
    colors = ['#F44336' if w<0 else '#4CAF50' for w in tbl['woe']]
    ax.bar(range(len(tbl)), tbl['woe'], color=colors, edgecolor='white')
    ax.axhline(y=0, color='black', lw=0.5)
    ax.set_title(f'{var} (IV={iv_results[var]:.3f})', fontweight='bold')
plt.tight_layout(); plt.show()

# --- WOE ENCODING COMPLETO ---
def apply_woe(df, feature, target, n_bins=10):
    """Aplica WOE encoding y retorna columna transformada."""
    temp = df[[feature, target]].dropna().copy()
    temp['bin'] = pd.qcut(temp[feature], q=n_bins, duplicates='drop')
    g = temp.groupby('bin', observed=True)[target].agg(['sum','count'])
    g.columns = ['n_d','tot']; g['n_nd'] = g['tot']-g['n_d']
    eps = 0.0001
    g['woe'] = np.log((g['n_nd']/g['n_nd'].sum()).clip(lower=eps) /
                      (g['n_d']/g['n_d'].sum()).clip(lower=eps))
    temp['woe_val'] = temp['bin'].map(dict(g['woe']))
    return temp['woe_val']

data_woe = data[['default']].copy()
for var in vars_std:
    data_woe[f'{var}_woe'] = apply_woe(data, var, 'default').values

# WOE categórica
g = data.groupby('tipo_empleo')['default'].agg(['sum','count'])
g.columns = ['n_d','tot']; g['n_nd'] = g['tot']-g['n_d']
eps = 0.0001
g['woe'] = np.log((g['n_nd']/g['n_nd'].sum()).clip(lower=eps) /
                  (g['n_d']/g['n_d'].sum()).clip(lower=eps))
data_woe['tipo_empleo_woe'] = data['tipo_empleo'].map(dict(g['woe']))

print("\nDataset WOE:")
print(data_woe.head(5).round(3))

# =============================================================================
# MÓDULO 6: SCORECARD COMPLETO
# Autor: Victor Miranda
# =============================================================================

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score, roc_curve, classification_report

# --- PASO 1: TRAIN/TEST SPLIT ---
woe_feats = [c for c in data_woe.columns if '_woe' in c]
X = data_woe[woe_feats].values
y = data_woe['default'].values

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=42, stratify=y)

print(f"Train: {X_train.shape[0]} | Test: {X_test.shape[0]}")
print(f"Default rate — Train: {y_train.mean():.2%}, Test: {y_test.mean():.2%}")

# --- PASO 2: REGRESIÓN LOGÍSTICA ---
model = LogisticRegression(C=1.0, solver='lbfgs', max_iter=1000, random_state=42)
model.fit(X_train, y_train)

coef_df = pd.DataFrame({'Variable': woe_feats,
                         'Coef': model.coef_[0]}).sort_values('Coef', ascending=False)
print(f"\nIntercepto: {model.intercept_[0]:.4f}")
print(coef_df.to_string(index=False))

# Verificar signos
neg = coef_df[coef_df['Coef'] < 0]
print(f"\n{'⚠ '+str(len(neg))+' coefs negativos' if len(neg)>0 else '✅ Todos positivos'}")

# --- PASO 3: EVALUACIÓN ---
y_prob_train = model.predict_proba(X_train)[:,1]
y_prob_test = model.predict_proba(X_test)[:,1]

auc_tr = roc_auc_score(y_train, y_prob_train)
auc_te = roc_auc_score(y_test, y_prob_test)
gini_tr, gini_te = 2*auc_tr-1, 2*auc_te-1

print(f"\nAUC  — Train: {auc_tr:.4f} | Test: {auc_te:.4f}")
print(f"Gini — Train: {gini_tr:.4f} | Test: {gini_te:.4f}")

# Curva ROC
fpr, tpr, thresholds = roc_curve(y_test, y_prob_test)

fig, axes = plt.subplots(1,2, figsize=(14,6))
axes[0].plot(fpr, tpr, 'b-', lw=2,
             label=f'AUC={auc_te:.3f} (Gini={gini_te:.3f})')
axes[0].plot([0,1],[0,1],'r--')
axes[0].set_xlabel('FPR'); axes[0].set_ylabel('TPR')
axes[0].set_title('Curva ROC — Test', fontweight='bold')
axes[0].legend()

# Youden optimal cutoff
youden = tpr - fpr
opt_idx = np.argmax(youden)
opt_thresh = thresholds[opt_idx]
print(f"\nCut-off (Youden): {opt_thresh:.4f}")
print(f"Sensitivity: {tpr[opt_idx]:.3f} | Specificity: {1-fpr[opt_idx]:.3f}")

# --- PASO 4: SCORECARD ---
pdo = 20; score_ref = 600; odds_ref = 50
factor = pdo / np.log(2)
offset = score_ref - factor * np.log(odds_ref)

print(f"\nSCORECARD: PDO={pdo}, Ref={score_ref}@{odds_ref}:1")
print(f"Factor={factor:.4f}, Offset={offset:.4f}")

# Calcular scores
log_odds = np.log((1-y_prob_test) / y_prob_test)
scores = np.round(offset + factor * log_odds).astype(int)

print(f"\nScores — Min: {scores.min()}, Max: {scores.max()}")
print(f"  Mean: {scores.mean():.0f}, Std: {scores.std():.0f}")

# Score distribution by group
axes[1].hist(scores[y_test==0], bins=40, alpha=0.6, color='#2196F3',
             density=True, label='No Default')
axes[1].hist(scores[y_test==1], bins=40, alpha=0.6, color='#F44336',
             density=True, label='Default')
axes[1].set_xlabel('Score'); axes[1].set_ylabel('Densidad')
axes[1].set_title('Distribución de Scores', fontweight='bold')
axes[1].legend()
plt.tight_layout(); plt.show()

# --- PASO 5: SCORECARD POR ATRIBUTO ---
n_vars = len(woe_feats)
intercept = model.intercept_[0]
coefs = model.coef_[0]

print("\nSCORECARD POR ATRIBUTO")
print("="*60)
for j, feat in enumerate(woe_feats):
    orig = feat.replace('_woe','')
    if orig in woe_tables:
        tbl = woe_tables[orig].copy()
        tbl['points'] = -(coefs[j] * tbl['woe'] * factor) + \
                       (offset - intercept * factor) / n_vars
        tbl['points'] = tbl['points'].round(0).astype(int)
        print(f"\n--- {orig} (β={coefs[j]:.4f}) ---")
        print(tbl[['bin','woe','def_rate','points']].to_string(index=False))

# --- PASO 6: ANÁLISIS POR BANDAS ---
score_df = pd.DataFrame({'score': scores, 'default': y_test})
score_df['band'] = pd.cut(score_df['score'], bins=10)

bands = score_df.groupby('band', observed=True).agg(
    n=('score','count'), n_def=('default','sum'),
    dr=('default','mean')).reset_index()
bands['pct'] = bands['n'] / bands['n'].sum()
bands['cum_dr'] = bands['n_def'].cumsum() / bands['n'].cumsum()

print("\nANÁLISIS POR BANDAS DE SCORE")
print(bands.to_string(index=False))