Driver turnover:
spotting the group about to tip - before the notices arrive

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

np.random.seed(303)

# --- Fahrer-Stammdaten ---
FAHRER_COUNT = 210
fahrer = []
for i in range(FAHRER_COUNT):
    alter = np.random.normal(42, 10)
    betriebszug = np.random.exponential(4.5)  # Jahre
    pendelzeit = np.random.choice([15,25,35,50,70],
                                   p=[0.2,0.3,0.25,0.15,0.1])
    schichttyp = np.random.choice(['Fernverkehr', 'Nahverkehr', 'Wechselbrücke'],
                                   p=[0.45, 0.35, 0.20])

    # Fluktuation-Wahrscheinlichkeit: Multi-Faktor-Modell
    base_churn = 0.12
    if betriebszug < 1.5: base_churn += 0.15       # Neue Fahrer gehen häufiger
    if pendelzeit > 50: base_churn += 0.10        # Langer Pendelweg
    if alter < 30: base_churn += 0.08             # Jüngere wechseln öfter
    if schichttyp == 'Fernverkehr': base_churn += 0.06

    hat_gekuendigt = np.random.random() < base_churn
    kuendigung_monat = np.random.randint(6, 24) if hat_gekuendigt else None

    fahrer.append({
        'fahrer_id': f'F-{i+1:03d}',
        'alter': int(np.clip(alter, 21, 64)),
        'betriebszugehoerigkeit_j': round(max(0.3, betriebszug), 1),
        'pendelzeit_min': pendelzeit,
        'schichttyp': schichttyp,
        'fuehrerschein_klasse': np.random.choice(['CE', 'C'], p=[0.85,0.15]),
        'hat_gekuendigt': hat_gekuendigt,
        'kuendigung_monat': kuendigung_monat,
    })

# --- Monatliche Betriebsdaten pro Fahrer (24 Monate) ---
records = []
for f in fahrer:
    for monat in range(1, 25):
        # Vorher-Signal: 3 Monate vor Kündigung ändern sich Muster
        vor_kuendigung = (f['hat_gekuendigt']
                          and f['kuendigung_monat'] - monat <= 3
                          and f['kuendigung_monat'] - monat >= 0)

        ueberstunden = np.random.normal(18, 8)
        kranktage = np.random.poisson(0.8)
        wochenend_schichten = np.random.poisson(2.5)
        mehrtagestouren = np.random.poisson(4) if f['schichttyp'] == 'Fernverkehr' else 0
        beschwerden = 0

        # --- Signals vor Kündigung ---
        if vor_kuendigung:
            kranktage += np.random.poisson(1.5)     # Mehr Krankmeldungen
            ueberstunden -= np.random.uniform(3,8)  # Weniger Bereitschaft
            beschwerden += np.random.choice([0,1,1,2])

        # Fahrverhalten aus Telematik
        harshbraking = np.random.poisson(12)
        geschw_verstoesse = np.random.poisson(3)
        if vor_kuendigung:
            harshbraking += np.random.poisson(5)   # Unachtsamkeit steigt
            geschw_verstoesse += np.random.poisson(2)

        is_target_month = (f['hat_gekuendigt']
                           and monat == f['kuendigung_monat'])

        records.append({
            'fahrer_id': f['fahrer_id'],
            'monat': monat,
            'alter': f['alter'],
            'betriebszug_j': f['betriebszugehoerigkeit_j'],
            'pendelzeit_min': f['pendelzeit_min'],
            'schichttyp': f['schichttyp'],
            'ueberstunden': round(max(0, ueberstunden), 1),
            'kranktage': kranktage,
            'wochenend_schichten': wochenend_schichten,
            'mehrtagestouren': mehrtagestouren,
            'beschwerden': beschwerden,
            'harshbraking_events': harshbraking,
            'geschw_verstoesse': geschw_verstoesse,
            'kuendigung': 1 if is_target_month else 0,
        })

df = pd.DataFrame(records)
# Beobachtung ueber 24 Monate -> Jahresrate
churn_jahr = df.groupby('fahrer_id')['kuendigung'].max().mean() / 2
print(f"Dataset: {len(df):,} Monatsdatensätze, {FAHRER_COUNT} Fahrer")
print(f"Fluktuation: {churn_jahr*100:.1f}% pro Jahr")

# Vergleich: Kündiger vs. Bleiber (3 Monate vor Austritt)
kuendiger = df[df['kuendigung'] == 1]['fahrer_id'].unique()

# Letzte 3 Monate vor Kündigung vs. normaler Betrieb
pre_churn = df[(df['fahrer_id'].isin(kuendiger))]
stable = df[~df['fahrer_id'].isin(kuendiger)]

compare = pd.DataFrame({
    'Kündiger (3M vorher)': pre_churn.groupby('fahrer_id').tail(3).mean(numeric_only=True),
    'Bleiber (Ø)': stable.mean(numeric_only=True),
})
print(compare.T[['ueberstunden','kranktage','wochenend_schichten',
                 'beschwerden','harshbraking_events']].round(1))

# Churn-Rate nach Betriebszugehörigkeit
fahrer_df = pd.DataFrame(fahrer)
bins = [0, 1, 2, 4, 8, 20]
fahrer_df['bz_gruppe'] = pd.cut(fahrer_df['betriebszugehoerigkeit_j'], bins)
churn_by_bz = fahrer_df.groupby('bz_gruppe')['hat_gekuendigt'].mean() * 100
print("\nChurn-Rate nach Betriebszugehörigkeit:")
print(churn_by_bz.round(1))

# --- Trend-Features: Veränderung über die letzten 3 Monate ---
for col in ['kranktage', 'ueberstunden', 'harshbraking_events', 'beschwerden']:
    df[f'{col}_trend_3m'] = df.groupby('fahrer_id')[col].transform(
        lambda x: x.rolling(3, min_periods=1).mean().diff()
    )

# --- Verhältnis-Features ---
df['ue_zu_we_ratio'] = df['ueberstunden'] / (df['wochenend_schichten'] + 1)
df['krank_zu_betrieb_ratio'] = df['kranktage'] / (df['betriebszug_j'] + 0.5)

# --- Volatilitäts-Feature: Wie unregelmäßig sind die Muster? ---
df['ueberstunden_volatilitaet'] = df.groupby('fahrer_id')['ueberstunden'].transform(
    lambda x: x.rolling(6, min_periods=2).std()
)

# --- Encoding ---
from sklearn.preprocessing import LabelEncoder
le = LabelEncoder()
df['schicht_enc'] = le.fit_transform(df['schichttyp'])

features = [
    # Stammdaten
    'alter', 'betriebszug_j', 'pendelzeit_min', 'schicht_enc',
    # Monatliche Betriebsdaten
    'ueberstunden', 'kranktage', 'wochenend_schichten',
    'mehrtagestouren', 'beschwerden',
    # Telematik
    'harshbraking_events', 'geschw_verstoesse',
    # Trend-Features
    'kranktage_trend_3m', 'ueberstunden_trend_3m',
    'harshbraking_events_trend_3m',
    # Verhältnis-Features
    'krank_zu_betrieb_ratio', 'ueberstunden_volatilitaet',
]

print(f"Feature-Matrix: {len(features)} Features × {len(df)} Monate")

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import TimeSeriesSplit
from sklearn.metrics import accuracy_score, recall_score
from sklearn.calibration import CalibratedClassifierCV

# --- Train/Test: Zeitlich (Monat 1-18 = Train, 19-24 = Test) ---
train = df[df['monat'] <= 18].dropna(subset=features)
test = df[df['monat'] > 18].dropna(subset=features)

# Target: Kündigt der Fahrer in den nächsten 3 Monaten?
y_train = train['kuendigung']
y_test = test['kuendigung']
X_train = train[features]
X_test = test[features]

# --- Random Forest mit Kalibrierung ---
rf = RandomForestClassifier(
    n_estimators=400,
    max_depth=8,
    min_samples_leaf=10,
    class_weight='balanced_subsample',
    random_state=42,
    n_jobs=-1
)
model = CalibratedClassifierCV(rf, cv=5, method='isotonic')
model.fit(X_train, y_train)

# --- Evaluation (Größenordnung Correll et al. 2024) ---
y_pred = model.predict(X_test)

print(f"Accuracy: {accuracy_score(y_test, y_pred):.2f}")
print(f"Recall:   {recall_score(y_test, y_pred):.2f}")

# Auswertung: Aggregation je Schichtgruppe - keine Einzel-Scores im Report

from lifelines import KaplanMeierFitter

# Survival-Kurve nach Schichttyp
kmf = KaplanMeierFitter()

fig, ax = plt.subplots(figsize=(10, 6))
for schicht in ['Fernverkehr', 'Nahverkehr', 'Wechselbrücke']:
    mask = fahrer_df['schichttyp'] == schicht
    durations = fahrer_df[mask]['betriebszugehoerigkeit_j'] * 12  # in Monaten
    events = fahrer_df[mask]['hat_gekuendigt'].astype(int)
    kmf.fit(durations, events, label=schicht)
    kmf.plot_survival_function(ax=ax)

ax.set_xlabel('Monate seit Eintritt')
ax.set_ylabel('Verbleibswahrscheinlichkeit')
ax.set_title('Kaplan-Meier: Überlebenskurve nach Schichttyp')

# Median-Überlebenszeit
for schicht in ['Fernverkehr', 'Nahverkehr', 'Wechselbrücke']:
    mask = fahrer_df['schichttyp'] == schicht
    churn_rate = fahrer_df[mask]['hat_gekuendigt'].mean()
    print(f"{schicht}: {churn_rate*100:.1f}% Fluktuation")

# === KOSTEN EINES FAHRERVERLUSTS (US-Band ~$6.000-12.000) ===
verlust_kosten = {
    'recruiting': 2500,         # Stellenanzeige, Agentur, Vorstellungsgespräche
    'einarbeitung': 2800,       # reduzierte Produktivität + Mentor
    'fahrzeug_leerstand': 2200, # LKW ohne Fahrer × Deckungsbeitrag
    'ueberstunden_team': 1200, # Kollegen übernehmen Touren
    'admin_hr': 800,            # Verwaltungsaufwand
}
kosten_pro_verlust = sum(verlust_kosten.values())  # €9.500

# === SZENARIO: Modell-gestützte Retention (Gruppenebene) ===
fahrer_gesamt = 210
fluktuation_jahr = 0.20  # Annahme - keine publizierte DE-Quote
abgaenge = int(fahrer_gesamt * fluktuation_jahr)  # 42

# Modell erkennt 60% vorab (Recall, Correll et al. 2024)
erkannt = int(abgaenge * 0.60)  # 25

# Interventionen halten 30% der erkannten Abgänge
gehalten = int(erkannt * 0.30)  # 7

ersparnis = gehalten * kosten_pro_verlust
kosten_retention = gehalten * 1800  # Ø Kosten pro Intervention
netto = ersparnis - kosten_retention

print(f"Abgänge/Jahr:                {abgaenge}")
print(f"Vom Modell vorab erkannt:    {erkannt}")
print(f"Durch Intervention gehalten: {gehalten}")
print(f"Kosten pro Verlust:          €{kosten_pro_verlust:,}")
print(f"Brutto-Ersparnis:            €{ersparnis:,}")
print(f"Interventionskosten:         €{kosten_retention:,}")
print(f"NETTO-ERSPARNIS:             €{netto:,}")