Building Factor Models on India Stack Data: UPI, AA, and GST as Alpha Signals

import numpy as np
import pandas as pd
from scipy import stats
from scipy.stats import norm
from typing import Dict, Tuple, Optional
import warnings
warnings.filterwarnings('ignore')


# ---------------------------------------------------------------------------
# 1. Data loading and preprocessing
# ---------------------------------------------------------------------------

def load_upi_data(filepath: str) -> pd.DataFrame:
    """
    Load NPCI UPI aggregate data. Expected columns:
    date, merchant_category, transaction_count, transaction_value, state
    """
    df = pd.read_csv(filepath, parse_dates=['date'])
    sector_map = {
        'grocery': 'FMCG',
        'fuel': 'Energy',
        'telecom_recharge': 'Telecom',
        'education': 'Education',
        'insurance': 'Insurance',
        'mutual_fund': 'Financial Services',
        'utilities': 'Utilities',
        'ecommerce': 'Consumer Discretionary',
        'travel': 'Travel & Hospitality',
        'healthcare': 'Pharma & Healthcare',
    }
    df['sector'] = df['merchant_category'].map(sector_map)
    df = df.dropna(subset=['sector'])
    return df.groupby(['date', 'sector'])['transaction_value'].sum().reset_index()


def load_gst_data(filepath: str) -> pd.DataFrame:
    """
    Load GSTN aggregate data. Expected columns:
    date, sector, output_tax, input_tax_credit, num_filers
    """
    df = pd.read_csv(filepath, parse_dates=['date'])
    df['itc_ratio'] = df['input_tax_credit'] / df['output_tax'].replace(0, np.nan)
    return df


def load_nse500_universe(filepath: str) -> pd.DataFrame:
    """
    Point-in-time NSE 500 membership with sector classifications.
    Columns: date, symbol, sector, market_cap, free_float_mcap, revenue_by_sector (JSON)
    """
    df = pd.read_csv(filepath, parse_dates=['date'])
    return df


# ---------------------------------------------------------------------------
# 2. Factor construction
# ---------------------------------------------------------------------------

def winsorize(series: pd.Series, limits: Tuple[float, float] = (0.01, 0.99)) -> pd.Series:
    """Winsorize at given percentiles to handle outliers."""
    lower = series.quantile(limits[0])
    upper = series.quantile(limits[1])
    return series.clip(lower, upper)


def cross_sectional_normalize(scores: pd.Series) -> pd.Series:
    """Rank-based inverse normal transformation."""
    n = len(scores)
    ranks = scores.rank(method='average')
    uniform = ranks / (n + 1)
    return pd.Series(norm.ppf(uniform), index=scores.index)


def construct_upi_factor(
    upi_data: pd.DataFrame,
    universe: pd.DataFrame,
    lookback: int = 12
) -> pd.DataFrame:
    """
    Construct UPI Consumption Velocity factor.
    Returns DataFrame with columns: date, symbol, upi_factor
    """
    # Compute sector-level z-scores
    sector_monthly = upi_data.pivot_table(
        index='date', columns='sector', values='transaction_value'
    )
    sector_monthly = sector_monthly.sort_index()

    rolling_mean = sector_monthly.rolling(lookback).mean()
    rolling_std = sector_monthly.rolling(lookback).std()
    sector_zscore = (sector_monthly - rolling_mean) / rolling_std

    # Map to stocks
    records = []
    for date in sector_zscore.index:
        if date not in universe['date'].values:
            continue
        month_universe = universe[universe['date'] == date]
        month_scores = sector_zscore.loc[date]

        for _, stock in month_universe.iterrows():
            sector = stock['sector']
            if sector in month_scores.index and not np.isnan(month_scores[sector]):
                records.append({
                    'date': date,
                    'symbol': stock['symbol'],
                    'upi_factor_raw': month_scores[sector]
                })

    result = pd.DataFrame(records)

    # Cross-sectional normalization per date
    result['upi_factor'] = result.groupby('date')['upi_factor_raw'].transform(
        lambda x: cross_sectional_normalize(winsorize(x))
    )
    return result[['date', 'symbol', 'upi_factor']]


def construct_gst_factor(
    gst_data: pd.DataFrame,
    universe: pd.DataFrame,
    ma_window: int = 6
) -> pd.DataFrame:
    """
    Construct GST Revenue Quality factor.
    Returns DataFrame with columns: date, symbol, gst_factor
    """
    sector_itc = gst_data.pivot_table(
        index='date', columns='sector', values='itc_ratio'
    ).sort_index()

    # Negative of deviation from MA: improving margins = positive score
    ma = sector_itc.rolling(ma_window).mean()
    gst_signal = -(sector_itc - ma)

    records = []
    for date in gst_signal.index:
        if date not in universe['date'].values:
            continue
        month_universe = universe[universe['date'] == date]
        month_scores = gst_signal.loc[date]

        for _, stock in month_universe.iterrows():
            sector = stock['sector']
            if sector in month_scores.index and not np.isnan(month_scores[sector]):
                records.append({
                    'date': date,
                    'symbol': stock['symbol'],
                    'gst_factor_raw': month_scores[sector]
                })

    result = pd.DataFrame(records)
    result['gst_factor'] = result.groupby('date')['gst_factor_raw'].transform(
        lambda x: cross_sectional_normalize(winsorize(x))
    )
    return result[['date', 'symbol', 'gst_factor']]


def construct_aa_credit_factor(
    credit_data: pd.DataFrame,
    delinquency_data: pd.DataFrame,
    universe: pd.DataFrame,
    penalty_lambda: float = 2.0
) -> pd.DataFrame:
    """
    Construct AA Credit Flow factor with delinquency adjustment.
    credit_data columns: date, sector, credit_outstanding
    delinquency_data columns: date, sector, delinquency_rate_90plus
    """
    credit_pivot = credit_data.pivot_table(
        index='date', columns='sector', values='credit_outstanding'
    ).sort_index()

    delinq_pivot = delinquency_data.pivot_table(
        index='date', columns='sector', values='delinquency_rate_90plus'
    ).sort_index()

    # Credit acceleration: YoY change in growth rate
    log_credit = np.log(credit_pivot)
    credit_growth = log_credit.diff(1)
    credit_accel = credit_growth - credit_growth.shift(12)

    # Delinquency ratio (current / year-ago)
    delinq_ratio = delinq_pivot / delinq_pivot.shift(12)

    # Quality-adjusted credit signal
    credit_signal = credit_accel * (1 - penalty_lambda * (delinq_ratio - 1).clip(lower=0))

    records = []
    for date in credit_signal.index:
        if date not in universe['date'].values:
            continue
        month_universe = universe[universe['date'] == date]
        month_scores = credit_signal.loc[date]

        for _, stock in month_universe.iterrows():
            sector = stock['sector']
            if sector in month_scores.index and not np.isnan(month_scores[sector]):
                records.append({
                    'date': date,
                    'symbol': stock['symbol'],
                    'aa_factor_raw': month_scores[sector]
                })

    result = pd.DataFrame(records)
    result['aa_factor'] = result.groupby('date')['aa_factor_raw'].transform(
        lambda x: cross_sectional_normalize(winsorize(x))
    )
    return result[['date', 'symbol', 'aa_factor']]


# ---------------------------------------------------------------------------
# 3. Orthogonalization pipeline
# ---------------------------------------------------------------------------

def orthogonalize_factor(
    novel_factor: pd.Series,
    control_factors: pd.DataFrame
) -> pd.Series:
    """
    Regress novel factor on control factors (Fama-French + momentum),
    return residuals as orthogonalized factor.
    """
    valid = novel_factor.notna() & control_factors.notna().all(axis=1)
    y = novel_factor[valid].values
    X = control_factors[valid].values
    X = np.column_stack([np.ones(len(X)), X])

    # OLS
    beta = np.linalg.lstsq(X, y, rcond=None)[0]
    residuals = y - X @ beta

    result = pd.Series(np.nan, index=novel_factor.index)
    result[valid] = residuals
    return result


def run_orthogonalization(
    merged: pd.DataFrame,
    novel_cols: list,
    control_cols: list = ['mkt_excess', 'size', 'value', 'momentum', 'profitability']
) -> pd.DataFrame:
    """
    Cross-sectional orthogonalization for each date.
    """
    orthogonalized = {col: [] for col in novel_cols}
    dates = merged['date'].unique()

    for date in dates:
        mask = merged['date'] == date
        month_data = merged[mask]

        if len(month_data) < 50:
            continue

        controls = month_data[control_cols]

        for col in novel_cols:
            resid = orthogonalize_factor(month_data[col], controls)
            orthogonalized[col].append(
                pd.DataFrame({'date': date, 'symbol': month_data['symbol'].values,
                              col + '_orth': resid.values})
            )

    result_frames = {}
    for col in novel_cols:
        result_frames[col] = pd.concat(orthogonalized[col], ignore_index=True)

    # Merge all orthogonalized factors
    base = result_frames[novel_cols[0]]
    for col in novel_cols[1:]:
        base = base.merge(result_frames[col], on=['date', 'symbol'], how='outer')

    return base


# ---------------------------------------------------------------------------
# 4. IC analysis and factor evaluation
# ---------------------------------------------------------------------------

def compute_ic_series(
    factor_scores: pd.DataFrame,
    returns: pd.DataFrame,
    factor_col: str,
    return_col: str = 'fwd_return_1m'
) -> pd.DataFrame:
    """
    Compute monthly Information Coefficient (Spearman rank correlation
    between factor scores and forward returns).
    """
    merged = factor_scores.merge(returns, on=['date', 'symbol'], how='inner')
    ic_by_date = merged.groupby('date').apply(
        lambda df: stats.spearmanr(df[factor_col], df[return_col])[0]
        if len(df) > 30 else np.nan
    ).dropna()

    ic_df = pd.DataFrame({
        'date': ic_by_date.index,
        'ic': ic_by_date.values
    })
    return ic_df


def ic_summary(ic_series: pd.DataFrame) -> Dict:
    """Summary statistics for IC series."""
    ic = ic_series['ic']
    return {
        'mean_ic': ic.mean(),
        'std_ic': ic.std(),
        'ic_ir': ic.mean() / ic.std() if ic.std() > 0 else 0,
        'pct_positive': (ic > 0).mean(),
        'hit_rate': (ic > 0).mean(),
        't_stat': ic.mean() / (ic.std() / np.sqrt(len(ic))) if ic.std() > 0 else 0,
    }


def ic_decay_analysis(
    factor_scores: pd.DataFrame,
    returns_multi: pd.DataFrame,
    factor_col: str,
    horizons: list = [1, 2, 3, 6, 12]
) -> pd.DataFrame:
    """
    Compute IC at multiple forward-return horizons to assess signal decay.
    returns_multi should have columns: fwd_return_1m, fwd_return_2m, etc.
    """
    results = []
    for h in horizons:
        ret_col = f'fwd_return_{h}m'
        if ret_col not in returns_multi.columns:
            continue
        ic_s = compute_ic_series(factor_scores, returns_multi, factor_col, ret_col)
        summary = ic_summary(ic_s)
        summary['horizon_months'] = h
        results.append(summary)
    return pd.DataFrame(results)


# ---------------------------------------------------------------------------
# 5. Fama-MacBeth regression for factor premium estimation
# ---------------------------------------------------------------------------

def fama_macbeth_regression(
    data: pd.DataFrame,
    return_col: str,
    factor_cols: list
) -> Dict:
    """
    Two-pass Fama-MacBeth regression.
    First pass: cross-sectional regression each period.
    Second pass: time-series average of cross-sectional coefficients.
    """
    dates = data['date'].unique()
    all_coeffs = []

    for date in sorted(dates):
        month = data[data['date'] == date].dropna(subset=[return_col] + factor_cols)
        if len(month) < 50:
            continue

        y = month[return_col].values
        X = month[factor_cols].values
        X = np.column_stack([np.ones(len(X)), X])

        try:
            beta = np.linalg.lstsq(X, y, rcond=None)[0]
            all_coeffs.append(beta)
        except np.linalg.LinAlgError:
            continue

    if not all_coeffs:
        return {}

    coeffs = np.array(all_coeffs)
    mean_coeffs = coeffs.mean(axis=0)
    std_coeffs = coeffs.std(axis=0)
    t_stats = mean_coeffs / (std_coeffs / np.sqrt(len(coeffs)))

    result = {'intercept': {'mean': mean_coeffs[0], 't_stat': t_stats[0]}}
    for i, col in enumerate(factor_cols):
        result[col] = {
            'mean_premium': mean_coeffs[i + 1],
            'std': std_coeffs[i + 1],
            't_stat': t_stats[i + 1],
            'significant_5pct': abs(t_stats[i + 1]) > 1.96
        }
    return result


# ---------------------------------------------------------------------------
# 6. Backtest engine
# ---------------------------------------------------------------------------

def construct_ls_portfolio(
    factor_scores: pd.DataFrame,
    factor_col: str,
    n_quantiles: int = 5,
    weight_col: Optional[str] = 'free_float_mcap'
) -> pd.DataFrame:
    """
    Construct long-short portfolio from factor quintiles.
    Long top quintile, short bottom quintile.
    """
    def assign_quantile(group):
        group['quantile'] = pd.qcut(
            group[factor_col], n_quantiles, labels=False, duplicates='drop'
        ) + 1
        return group

    scored = factor_scores.groupby('date', group_keys=False).apply(assign_quantile)

    # Long top quintile, short bottom quintile
    long_mask = scored['quantile'] == n_quantiles
    short_mask = scored['quantile'] == 1

    scored['ls_weight'] = 0.0

    if weight_col and weight_col in scored.columns:
        for date in scored['date'].unique():
            date_mask = scored['date'] == date
            long_stocks = date_mask & long_mask
            short_stocks = date_mask & short_mask

            long_weights = scored.loc[long_stocks, weight_col]
            short_weights = scored.loc[short_stocks, weight_col]

            if long_weights.sum() > 0:
                scored.loc[long_stocks, 'ls_weight'] = (
                    long_weights / long_weights.sum()
                )
            if short_weights.sum() > 0:
                scored.loc[short_stocks, 'ls_weight'] = (
                    -short_weights / short_weights.sum()
                )
    else:
        for date in scored['date'].unique():
            date_mask = scored['date'] == date
            n_long = (date_mask & long_mask).sum()
            n_short = (date_mask & short_mask).sum()

            if n_long > 0:
                scored.loc[date_mask & long_mask, 'ls_weight'] = 1.0 / n_long
            if n_short > 0:
                scored.loc[date_mask & short_mask, 'ls_weight'] = -1.0 / n_short

    return scored


def backtest_factor(
    portfolio: pd.DataFrame,
    returns: pd.DataFrame,
    return_col: str = 'fwd_return_1m',
    transaction_cost_bps: float = 30.0
) -> pd.DataFrame:
    """
    Backtest long-short factor portfolio with transaction costs.
    """
    merged = portfolio.merge(returns, on=['date', 'symbol'], how='inner')
    merged['weighted_return'] = merged['ls_weight'] * merged[return_col]

    # Portfolio return per period
    port_returns = merged.groupby('date')['weighted_return'].sum().reset_index()
    port_returns.columns = ['date', 'gross_return']

    # Estimate turnover from weight changes
    weights_wide = portfolio.pivot_table(
        index='date', columns='symbol', values='ls_weight', fill_value=0
    )
    turnover = weights_wide.diff().abs().sum(axis=1) / 2
    turnover_df = turnover.reset_index()
    turnover_df.columns = ['date', 'turnover']

    port_returns = port_returns.merge(turnover_df, on='date', how='left')
    port_returns['turnover'] = port_returns['turnover'].fillna(0)

    # Net return after transaction costs
    port_returns['tc'] = port_returns['turnover'] * transaction_cost_bps / 10000
    port_returns['net_return'] = port_returns['gross_return'] - port_returns['tc']

    # Cumulative returns
    port_returns['cum_gross'] = (1 + port_returns['gross_return']).cumprod()
    port_returns['cum_net'] = (1 + port_returns['net_return']).cumprod()

    return port_returns


def backtest_summary(port_returns: pd.DataFrame, periods_per_year: int = 12) -> Dict:
    """Compute backtest performance statistics."""
    net = port_returns['net_return']
    cum = port_returns['cum_net']

    ann_return = (cum.iloc[-1]) ** (periods_per_year / len(net)) - 1
    ann_vol = net.std() * np.sqrt(periods_per_year)
    sharpe = ann_return / ann_vol if ann_vol > 0 else 0

    # Max drawdown
    rolling_max = cum.cummax()
    drawdown = (cum - rolling_max) / rolling_max
    max_dd = drawdown.min()

    return {
        'annualized_return': f'{ann_return:.2%}',
        'annualized_volatility': f'{ann_vol:.2%}',
        'sharpe_ratio': round(sharpe, 2),
        'max_drawdown': f'{max_dd:.2%}',
        'mean_monthly_turnover': f'{port_returns["turnover"].mean():.2%}',
        'total_months': len(net),
        'pct_positive_months': f'{(net > 0).mean():.1%}',
    }


# ---------------------------------------------------------------------------
# 7. Full pipeline
# ---------------------------------------------------------------------------

def run_full_pipeline(config: Dict) -> Dict:
    """
    End-to-end pipeline: load data, construct factors, orthogonalize,
    run Fama-MacBeth, backtest, and produce IC analysis.
    """
    # Load data
    upi = load_upi_data(config['upi_path'])
    gst = load_gst_data(config['gst_path'])
    universe = load_nse500_universe(config['universe_path'])
    returns = pd.read_csv(config['returns_path'], parse_dates=['date'])

    # Construct factors
    upi_factor = construct_upi_factor(upi, universe)
    gst_factor = construct_gst_factor(gst, universe)

    # Merge all factors with returns
    merged = universe[['date', 'symbol', 'free_float_mcap']].copy()
    merged = merged.merge(upi_factor, on=['date', 'symbol'], how='left')
    merged = merged.merge(gst_factor, on=['date', 'symbol'], how='left')
    merged = merged.merge(returns, on=['date', 'symbol'], how='left')

    # Orthogonalize against standard factors
    if all(c in merged.columns for c in ['mkt_excess', 'size', 'value', 'momentum']):
        orth = run_orthogonalization(
            merged, ['upi_factor', 'gst_factor'],
            ['mkt_excess', 'size', 'value', 'momentum']
        )
        merged = merged.merge(orth, on=['date', 'symbol'], how='left')

    # Fama-MacBeth
    novel_cols = ['upi_factor', 'gst_factor']
    available_cols = [c for c in novel_cols if c in merged.columns]
    fm_results = fama_macbeth_regression(merged, 'fwd_return_1m', available_cols)

    # IC analysis
    ic_results = {}
    for col in available_cols:
        ic_s = compute_ic_series(merged, returns, col)
        ic_results[col] = ic_summary(ic_s)

    # Backtest each factor
    bt_results = {}
    for col in available_cols:
        portfolio = construct_ls_portfolio(merged.dropna(subset=[col]), col)
        bt = backtest_factor(portfolio, returns)
        bt_results[col] = backtest_summary(bt)

    return {
        'fama_macbeth': fm_results,
        'ic_analysis': ic_results,
        'backtest': bt_results,
    }