""" Supply Chain Payment & Fulfillment Scoring — MARS Implementation Brandon Wilhoite Models payment likelihood (0-200 scale) and fulfillment likelihood (0-200 scale) at the invoice/PO level using Multivariate Adaptive Regression Splines (MARS) via py-earth. Score interpretation: 200 = paid/fulfilled as early as possible 100 = paid/fulfilled exactly on time 0 = 180+ days late or complete failure to pay/fulfill Tiers: 1. Invoice-level payment score (Company as Buyer) 2. PO-level fulfillment score (Company as Seller) 3. Company Buyer Profile (aggregated from invoice-level) 4. Company Seller Profile (aggregated from PO-level) 5. Company Overall Score (average of Buyer + Seller profiles) 6. Normalized scores (geo / industry / size peer-relative) Architecture note: All tiers share the same variable universe. MARS pruning selects relevant variables and weights per tier based on the outcome being predicted. Missingness is encoded as a binary feature (distinct from a score of zero). No imputation — absent data is semantically different from zero. """ import numpy as np import pandas as pd from pyearth import Earth from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split from sklearn.metrics import mean_squared_error, r2_score from scipy.stats import pearsonr import warnings warnings.filterwarnings('ignore') # ── CONSTANTS ───────────────────────────────────────────────────────────────── SCORE_MIN = 0 SCORE_MAX = 200 SCORE_ONTIME = 100 # paid exactly on due date SCORE_FLOOR = 0 # 180+ days late or no payment SCORE_CEILING = 200 # paid as early as contractually possible # Normalization cohort dimensions GEO_COHORTS = ['North America', 'EU/UK', 'LATAM', 'Africa', 'Middle East', 'Asia', 'Australia'] SIZE_BANDS = ['micro', 'small', 'medium', 'large'] # by revenue or employee count # ── VARIABLE DEFINITIONS ────────────────────────────────────────────────────── # Full variable universe — same set feeds all tiers; MARS selects per outcome. # Variables marked [BINARY_IF_MISSING] get a companion missingness flag. FINANCIAL_RATIO_VARS = [ 'current_ratio', 'quick_ratio', 'cash_ratio', 'dscr', # Debt Service Coverage Ratio 'debt_service_coverage_ratio', 'fixed_charge_coverage_ratio', 'coverage_ratio', 'solvency_ratio', 'leverage', 'gross_profit_margin', 'working_capital', 'free_cash_flow', 'discretionary_cash_flow', 'net_operating_profit', 'operating_income', 'net_profit', 'ebitda', 'cash_conversion_cycle', 'receivable_turnover', ] BALANCE_SHEET_VARS = [ 'total_assets', 'current_liabilities', 'long_term_liabilities', 'liabilities', 'equity', 'shareholder_equity', 'retained_earnings', 'intangible_assets', 'tangible_asset_pct', 'asset_pct', 'asset_cr', # Asset Coverage Ratio 'cash', 'cash_on_hand', 'gross_sales', 'revenue', 'cogs', 'cost_of_sales', 'market_cap', 'pe_ratio', 'market_value_of_equity', ] PAYMENT_BEHAVIOR_VARS = [ 'avg_dbt', # Average Days Beyond Term 'high_dbt', # Highest recorded DBT 'days_beyond_term', # Current DBT 'prev_transaction_dbt', # DBT on most recent transaction 'prev_x_transactions_dbt', # Rolling average across recent transactions 'prev_90_days_dbt', 'prev_180_days_dbt', 'prev_year_dbt', 'avg_credit', 'high_credit', # Highest credit line extended 'credit_inquiries', 'num_delinquencies', 'num_bankruptcies', 'recurring_legal_cost', 'credit_revenue_ratio', # Credit extended relative to revenue 'total_spend', ] RISK_SCORE_VARS = [ 'probability_of_default', 'probability_of_delinquency', 'probability_of_failure', 'lgd', # Loss Given Default 'altman_z_score', 'internal_payscore', # Prior internal score if exists ] COMPANY_PROFILE_VARS = [ 'employee_count', 'years_established', 'revenue_growth', 'roa', 'roe', 'supplier_count', 'receivable_count', 'esg_flag', # ESG self-identification flag — negative predictor ] # Categorical variables — one-hot encoded before fitting CATEGORICAL_VARS = [ 'industry', 'region', 'geo', 'country', ] # Binary missingness flags — generated automatically for these variables # A company with no credit score is categorically different from a score of 0 BINARY_IF_MISSING = [ 'probability_of_default', 'probability_of_delinquency', 'probability_of_failure', 'altman_z_score', 'avg_dbt', 'high_credit', 'market_cap', 'ebitda', 'revenue', ] ALL_NUMERIC_VARS = ( FINANCIAL_RATIO_VARS + BALANCE_SHEET_VARS + PAYMENT_BEHAVIOR_VARS + RISK_SCORE_VARS + COMPANY_PROFILE_VARS ) # ── FEATURE ENGINEERING ─────────────────────────────────────────────────────── def add_missingness_flags(df: pd.DataFrame) -> pd.DataFrame: """ For variables in BINARY_IF_MISSING, add a companion binary column indicating whether the value was originally missing. Absent bureau data ≠ a score of zero. A company with no credit history is a distinct risk category from a company with a history score of 0. The missingness flag lets MARS treat these as separate cases. """ for var in BINARY_IF_MISSING: if var in df.columns: flag_col = f'{var}_is_missing' df[flag_col] = df[var].isna().astype(int) return df def encode_categoricals(df: pd.DataFrame) -> pd.DataFrame: """One-hot encode categorical variables, drop first to avoid multicollinearity.""" for var in CATEGORICAL_VARS: if var in df.columns: dummies = pd.get_dummies(df[var], prefix=var, drop_first=True) df = pd.concat([df.drop(columns=[var]), dummies], axis=1) return df def encode_binary_flags(df: pd.DataFrame) -> pd.DataFrame: """ Convert boolean/string binary fields to 0/1. has_bankruptcies, has_delinquencies treated as hard negative signals. esg_flag: empirically a negative payment predictor — companies self-identifying ESG focus showed systematically worse payment timing. """ binary_cols = ['has_bankruptcies', 'has_delinquencies', 'esg_flag'] for col in binary_cols: if col in df.columns: df[col] = df[col].map({True: 1, False: 0, 'yes': 1, 'no': 0, 1: 1, 0: 0}) return df def impute_with_missingness_flag(df: pd.DataFrame) -> pd.DataFrame: """ After flagging missingness, fill NaN with column median. The median fill is neutral — the missingness flag carries the signal. """ numeric_cols = df.select_dtypes(include=[np.number]).columns for col in numeric_cols: if df[col].isna().any(): df[col] = df[col].fillna(df[col].median()) return df def clip_score(score: float) -> float: """Enforce [0, 200] bounds on output.""" return float(np.clip(score, SCORE_MIN, SCORE_MAX)) def prepare_features(df: pd.DataFrame) -> pd.DataFrame: """Full preprocessing pipeline.""" df = df.copy() df = add_missingness_flags(df) df = encode_binary_flags(df) df = encode_categoricals(df) df = impute_with_missingness_flag(df) return df # ── SCORE CONSTRUCTION ──────────────────────────────────────────────────────── def invoice_score_from_timing(days_beyond_term: float, paid: bool) -> float: """ Convert observed payment timing to [0, 200] score. Score interpretation: 200 = paid as early as contractually possible 100 = paid exactly on due date (DBT = 0) 0 = 180+ days late or no payment Early payment (negative DBT) scores above 100. The exact early-payment ceiling and late-payment floor logic reflects the design intent without reproducing patented specifics. Args: days_beyond_term: actual DBT observed (negative = early) paid: whether invoice was ever paid Returns: float score in [0, 200] """ if not paid: return SCORE_FLOOR if days_beyond_term <= 0: # Early or on-time payment — scale from 100 to 200 # Earlier payment = higher score, capped at 200 early_bonus = min(abs(days_beyond_term), 90) / 90 * 100 return clip_score(SCORE_ONTIME + early_bonus) else: # Late payment — scale from 100 down to 0 over 180 days late_penalty = min(days_beyond_term, 180) / 180 * 100 return clip_score(SCORE_ONTIME - late_penalty) # ── MARS MODEL ──────────────────────────────────────────────────────────────── class PaymentScoreModel: """ MARS model for predicting invoice-level payment likelihood. Predicts a continuous score [0, 200] representing payment timing quality. Uses py-earth's Earth implementation of MARS, which automatically: - Identifies non-linear relationships via hinge functions - Discovers interaction terms between variables - Prunes non-contributory terms via GCV (generalized cross-validation) Why MARS over logistic regression or tree methods: - Payment behavior has threshold effects (5 days late ≠ 30 days late, but 120 days late ≈ 150 days late in terms of risk signal) - Industry and geography create structural breaks that must be discovered from data, not pre-specified - MARS hinge functions are interpretable and auditable — important for regulatory and underwriter review - Unlike black-box methods, MARS can explain why a score changed """ def __init__(self, max_degree: int = 2, max_terms: int = 50, min_samples_leaf: int = 10): """ Args: max_degree: maximum interaction degree between variables. 2 = pairwise interactions (recommended starting point) max_terms: maximum basis functions before pruning. Higher = more complex model, higher overfitting risk. min_samples_leaf: minimum observations per MARS leaf node. Higher = more conservative, less overfitting. """ self.max_degree = max_degree self.max_terms = max_terms self.min_samples_leaf = min_samples_leaf self.model = None self.scaler = StandardScaler() self.feature_names = None self.is_fitted = False def fit(self, X: pd.DataFrame, y: pd.Series) -> 'PaymentScoreModel': """ Fit MARS model. Args: X: feature matrix, preprocessed (missingness flags added, categoricals encoded, no NaN) y: target scores [0, 200] """ self.feature_names = list(X.columns) X_scaled = self.scaler.fit_transform(X) self.model = Earth( max_degree=self.max_degree, max_terms=self.max_terms, min_samples_leaf=self.min_samples_leaf, allow_linear=True, # allow linear terms alongside hinges enable_pruning=True, # GCV-based backward pruning ) self.model.fit(X_scaled, y) self.is_fitted = True return self def predict(self, X: pd.DataFrame) -> np.ndarray: """Predict scores, clipped to [0, 200].""" if not self.is_fitted: raise RuntimeError("Model must be fitted before predicting.") X_scaled = self.scaler.transform(X[self.feature_names]) raw = self.model.predict(X_scaled) return np.array([clip_score(s) for s in raw]) def feature_importance(self) -> pd.DataFrame: """ Return variable importance from MARS model. MARS importance is based on the GCV improvement attributable to each variable across all basis functions it appears in. Higher = more predictive. """ if not self.is_fitted: raise RuntimeError("Model must be fitted first.") importances = self.model.feature_importances_ return ( pd.DataFrame({ 'variable': self.feature_names, 'importance': importances }) .sort_values('importance', ascending=False) .reset_index(drop=True) ) def summary(self) -> str: """Return py-earth model summary string.""" if not self.is_fitted: raise RuntimeError("Model must be fitted first.") return self.model.summary() class FulfillmentScoreModel(PaymentScoreModel): """ MARS model for predicting PO-level fulfillment likelihood. Identical architecture to PaymentScoreModel — same variable universe, same [0, 200] scale, same MARS implementation. Different outcome: fulfillment timing and completeness rather than payment timing. 200 = fulfilled ahead of schedule, complete and accurate 100 = fulfilled exactly on time, complete 0 = 180+ days late, unfulfilled, or significant chargeback/return MARS selects different variable weights for fulfillment vs. payment because the behavioral drivers differ: seller quality metrics (chargeback rates, return rates, order accuracy) dominate here while cash flow ratios dominate payment scoring. """ pass # Identical architecture, outcome variable differs at training time # ── TIER AGGREGATION ───────────────────────────────────────────────────────── def aggregate_to_company_profile( invoice_scores: pd.Series, company_ids: pd.Series, weights: pd.Series = None ) -> pd.Series: """ Aggregate invoice-level scores to company profile score. Dollar-weighted aggregation: larger invoices carry proportionally more weight in the company profile than small ones. Args: invoice_scores: predicted score per invoice [0, 200] company_ids: company identifier per invoice weights: invoice amounts for dollar-weighting (optional; equal weighting if None) Returns: Series indexed by company_id with aggregated profile score """ df = pd.DataFrame({ 'company_id': company_ids, 'score': invoice_scores, 'weight': weights if weights is not None else pd.Series( np.ones(len(invoice_scores)), index=invoice_scores.index ) }) def weighted_mean(group): return np.average(group['score'], weights=group['weight']) return df.groupby('company_id').apply(weighted_mean) def company_overall_score( buyer_profile: pd.Series, seller_profile: pd.Series ) -> pd.Series: """ Company overall score = strict average of buyer and seller profiles. Both dimensions are equally weighted by design: a company's payment reliability (buyer) is no more or less important than its fulfillment reliability (seller) in a supply chain context. Companies operating only as buyer or only as seller receive their single available profile score as the overall score. """ combined = pd.DataFrame({ 'buyer': buyer_profile, 'seller': seller_profile }) # Where only one profile exists, use that profile as the overall score return combined.mean(axis=1, skipna=True) # ── NORMALIZATION ───────────────────────────────────────────────────────────── def normalize_within_cohort( scores: pd.Series, cohort_labels: pd.Series, method: str = 'percentile' ) -> pd.Series: """ Normalize scores within a peer cohort (geography, industry, or size band). Peer-relative scoring addresses a fundamental bias in absolute credit scoring: a company with a score of 55 may be the strongest payer in its geography and industry, or it may be below average. Absolute scores obscure this context. Args: scores: raw company overall scores [0, 200] cohort_labels: cohort identifier per company (geo, industry, or size) method: 'percentile' (rank within cohort, scaled to [0, 200]) or 'zscore' (standardized within cohort) Returns: Normalized scores in [0, 200] """ normalized = scores.copy().astype(float) for cohort in cohort_labels.unique(): mask = cohort_labels == cohort cohort_scores = scores[mask] if len(cohort_scores) < 2: # Insufficient cohort size — retain absolute score continue if method == 'percentile': # Rank within cohort, scale to [0, 200] ranks = cohort_scores.rank(pct=True) normalized[mask] = ranks * SCORE_MAX elif method == 'zscore': # Z-score within cohort, then map to [0, 200] via sigmoid z = (cohort_scores - cohort_scores.mean()) / cohort_scores.std() # Sigmoid maps z-scores to (0, 1), scale to [0, 200] sigmoid = 1 / (1 + np.exp(-z)) normalized[mask] = sigmoid * SCORE_MAX return normalized.clip(SCORE_MIN, SCORE_MAX) def normalize_all_dimensions( overall_scores: pd.Series, geo_labels: pd.Series, industry_labels: pd.Series, size_labels: pd.Series, method: str = 'percentile' ) -> pd.DataFrame: """ Compute all normalization dimensions for the overall company score. Returns a DataFrame with columns: score_global — raw overall score score_geo — normalized within geographic region score_industry — normalized within industry cohort score_size — normalized within company size band The global score is the primary output. Normalized dimensions provide interpretive context — they answer "how does this company compare to its peers" rather than overriding the behavioral signal. """ return pd.DataFrame({ 'score_global': overall_scores.clip(SCORE_MIN, SCORE_MAX), 'score_geo': normalize_within_cohort(overall_scores, geo_labels, method), 'score_industry': normalize_within_cohort(overall_scores, industry_labels, method), 'score_size': normalize_within_cohort(overall_scores, size_labels, method), }) # ── EVALUATION ──────────────────────────────────────────────────────────────── def evaluate_model( y_true: np.ndarray, y_pred: np.ndarray, model_name: str = "Model" ) -> dict: """ Evaluate model performance against observed outcomes. Primary metric: Pearson correlation to actual payment behavior. This is the same metric used to evaluate bureau data quality — allowing direct comparison between bureau signals and the transactional model on the same scale. """ corr, p_value = pearsonr(y_true, y_pred) rmse = np.sqrt(mean_squared_error(y_true, y_pred)) r2 = r2_score(y_true, y_pred) results = { 'model': model_name, 'correlation': round(corr, 4), 'p_value': round(p_value, 6), 'rmse': round(rmse, 2), 'r2': round(r2, 4), 'n': len(y_true) } print(f"\n{'─' * 50}") print(f" {model_name}") print(f"{'─' * 50}") print(f" Correlation to actual payment behavior: {corr:.4f}") print(f" R²: {r2:.4f}") print(f" RMSE: {rmse:.2f}") print(f" p-value: {p_value:.6f}") print(f" n (invoices): {len(y_true)}") return results # ── PIPELINE ───────────────────────────────────────────────────────────────── def run_scoring_pipeline( invoice_df: pd.DataFrame, po_df: pd.DataFrame, company_df: pd.DataFrame, fit_payment_model: bool = True, fit_fulfillment_model: bool = True, test_size: float = 0.2, random_state: int = 42 ) -> dict: """ Full scoring pipeline from raw data to normalized company scores. Args: invoice_df: Invoice-level data. Must contain all variable columns plus 'company_id', 'invoice_amount', 'payment_score' (observed outcome, [0, 200]). po_df: PO-level data. Same structure with 'fulfillment_score' as outcome. company_df: Company-level metadata. Must contain 'company_id', 'geo', 'industry', 'size_band'. fit_payment_model: Whether to train payment model. fit_fulfillment_model: Whether to train fulfillment model. test_size: Proportion of data held out for evaluation. random_state: Reproducibility seed. Returns: dict with keys: 'payment_model': fitted PaymentScoreModel 'fulfillment_model': fitted FulfillmentScoreModel 'company_scores': DataFrame of normalized company scores 'evaluation': dict of evaluation metrics per tier """ results = {} # ── PAYMENT MODEL (Company as Buyer) ────────────────────────────────── if fit_payment_model and invoice_df is not None: print("\n[1/4] Training payment model (Company as Buyer)...") X_inv = prepare_features(invoice_df.drop( columns=['company_id', 'invoice_amount', 'payment_score'], errors='ignore' )) y_inv = invoice_df['payment_score'] X_train, X_test, y_train, y_test = train_test_split( X_inv, y_inv, test_size=test_size, random_state=random_state ) payment_model = PaymentScoreModel(max_degree=2, max_terms=50) payment_model.fit(X_train, y_train) results['payment_model'] = payment_model y_pred_inv = payment_model.predict(X_test) results['eval_payment'] = evaluate_model( y_test.values, y_pred_inv, "Payment Model (Invoice Level)" ) # Aggregate to buyer profile all_payment_scores = pd.Series( payment_model.predict(X_inv), index=invoice_df.index ) buyer_profiles = aggregate_to_company_profile( all_payment_scores, invoice_df['company_id'], weights=invoice_df.get('invoice_amount') ) print(f" Buyer profiles computed for {len(buyer_profiles)} companies.") else: buyer_profiles = pd.Series(dtype=float) # ── FULFILLMENT MODEL (Company as Seller) ───────────────────────────── if fit_fulfillment_model and po_df is not None: print("\n[2/4] Training fulfillment model (Company as Seller)...") X_po = prepare_features(po_df.drop( columns=['company_id', 'po_amount', 'fulfillment_score'], errors='ignore' )) y_po = po_df['fulfillment_score'] X_train_po, X_test_po, y_train_po, y_test_po = train_test_split( X_po, y_po, test_size=test_size, random_state=random_state ) fulfillment_model = FulfillmentScoreModel(max_degree=2, max_terms=50) fulfillment_model.fit(X_train_po, y_train_po) results['fulfillment_model'] = fulfillment_model y_pred_po = fulfillment_model.predict(X_test_po) results['eval_fulfillment'] = evaluate_model( y_test_po.values, y_pred_po, "Fulfillment Model (PO Level)" ) all_fulfillment_scores = pd.Series( fulfillment_model.predict(X_po), index=po_df.index ) seller_profiles = aggregate_to_company_profile( all_fulfillment_scores, po_df['company_id'], weights=po_df.get('po_amount') ) print(f" Seller profiles computed for {len(seller_profiles)} companies.") else: seller_profiles = pd.Series(dtype=float) # ── COMPANY OVERALL SCORE ───────────────────────────────────────────── print("\n[3/4] Computing company overall scores...") overall = company_overall_score(buyer_profiles, seller_profiles) print(f" Overall scores computed for {len(overall)} companies.") # ── NORMALIZATION ───────────────────────────────────────────────────── print("\n[4/4] Normalizing scores within cohorts...") # Align company metadata to scored companies meta = company_df.set_index('company_id').reindex(overall.index) normalized_scores = normalize_all_dimensions( overall_scores=overall, geo_labels=meta['geo'].fillna('Unknown'), industry_labels=meta['industry'].fillna('Unknown'), size_labels=meta['size_band'].fillna('Unknown'), method='percentile' ) # Attach buyer and seller profiles normalized_scores['score_buyer_profile'] = buyer_profiles.reindex(overall.index) normalized_scores['score_seller_profile'] = seller_profiles.reindex(overall.index) normalized_scores.index.name = 'company_id' results['company_scores'] = normalized_scores print(f"\n{'═' * 50}") print(f" Pipeline complete.") print(f" Companies scored: {len(normalized_scores)}") print(f" Score range (global): " f"{normalized_scores['score_global'].min():.1f} – " f"{normalized_scores['score_global'].max():.1f}") print(f"{'═' * 50}\n") return results # ── USAGE EXAMPLE ───────────────────────────────────────────────────────────── if __name__ == "__main__": """ Illustrative usage with synthetic data. Replace with actual invoice_df, po_df, company_df from your data source. """ np.random.seed(42) n_invoices = 1000 n_companies = 80 company_ids = np.random.choice([f"CO_{i:03d}" for i in range(n_companies)], n_invoices) # Synthetic invoice data — in production, sourced from Codat accounting API invoice_df = pd.DataFrame({ 'company_id': company_ids, 'invoice_amount': np.random.lognormal(8, 1.5, n_invoices), 'avg_dbt': np.random.normal(12, 20, n_invoices), 'high_dbt': np.random.normal(30, 40, n_invoices).clip(0), 'current_ratio': np.random.normal(1.5, 0.5, n_invoices).clip(0.1), 'dscr': np.random.normal(1.2, 0.4, n_invoices), 'revenue': np.random.lognormal(12, 1.5, n_invoices), 'ebitda': np.random.normal(50000, 80000, n_invoices), 'esg_flag': np.random.choice([0, 1], n_invoices, p=[0.85, 0.15]), 'industry': np.random.choice(['Manufacturing', 'Retail', 'Services', 'Tech'], n_invoices), 'geo': np.random.choice(GEO_COHORTS[:4], n_invoices), 'num_delinquencies': np.random.poisson(0.3, n_invoices), 'years_established': np.random.randint(1, 30, n_invoices), 'employee_count': np.random.lognormal(3, 1, n_invoices).astype(int), # Simulate some missing bureau data 'probability_of_default': np.where( np.random.random(n_invoices) > 0.15, np.random.beta(2, 5, n_invoices) * 100, np.nan ), }) # Synthetic payment scores — ESG companies score worse (empirical finding) base_score = ( 100 - invoice_df['avg_dbt'].clip(0, 90) * 0.5 - invoice_df['num_delinquencies'] * 8 + (invoice_df['current_ratio'] - 1).clip(-1, 1) * 15 - invoice_df['esg_flag'] * 25 # ESG flag = negative payment predictor + np.random.normal(0, 10, n_invoices) ).clip(SCORE_MIN, SCORE_MAX) invoice_df['payment_score'] = base_score # Synthetic PO data po_df = invoice_df.copy().rename( columns={'invoice_amount': 'po_amount', 'payment_score': 'fulfillment_score'} ) po_df['fulfillment_score'] = ( base_score + np.random.normal(0, 5, n_invoices) ).clip(SCORE_MIN, SCORE_MAX) # Company metadata company_df = pd.DataFrame({ 'company_id': [f"CO_{i:03d}" for i in range(n_companies)], 'geo': np.random.choice(GEO_COHORTS[:4], n_companies), 'industry': np.random.choice(['Manufacturing', 'Retail', 'Services', 'Tech'], n_companies), 'size_band': np.random.choice(SIZE_BANDS, n_companies), }) # Run pipeline results = run_scoring_pipeline( invoice_df=invoice_df, po_df=po_df, company_df=company_df ) # Show top/bottom companies by global score scores = results['company_scores'].sort_values('score_global', ascending=False) print("\nTop 5 companies by global score:") print(scores[['score_global', 'score_buyer_profile', 'score_seller_profile', 'score_geo']].head()) print("\nBottom 5 companies by global score:") print(scores[['score_global', 'score_buyer_profile', 'score_seller_profile', 'score_geo']].tail()) # Feature importance from payment model print("\nTop 10 variables by importance (Payment Model):") print(results['payment_model'].feature_importance().head(10))