fidc-backtest-engine/crates/fidc-core/src/metrics.rs

use std::collections::BTreeMap;

use chrono::{Datelike, NaiveDate};
use serde::{Deserialize, Serialize};

use crate::engine::DailyEquityPoint;
use crate::events::FillEvent;
use crate::portfolio::HoldingSummary;

const TRADING_DAYS_PER_YEAR: f64 = 252.0;
const MONTHS_PER_YEAR: f64 = 12.0;
const DEFAULT_RISK_FREE_RATE: f64 = 0.022;

#[derive(Debug, Clone, Default, Serialize, Deserialize)]
pub struct BacktestMetrics {
    pub total_return: f64,
    pub annual_return: f64,
    pub sharpe: f64,
    pub max_drawdown: f64,
    pub win_rate: f64,
    pub alpha: f64,
    pub beta: f64,
    pub benchmark_cumulative_return: f64,
    pub benchmark_net_value: f64,
    pub risk_free_rate: f64,
    pub monthly_excess_win_rate: f64,
    pub excess_cumulative_return: f64,
    pub excess_annual_return: f64,
    pub max_drawdown_duration_days: usize,
    pub total_trade_days: usize,
    pub sortino: f64,
    pub information_ratio: f64,
    pub tracking_error: f64,
    pub volatility: f64,
    pub excess_return: f64,
    pub excess_sharpe: f64,
    pub excess_volatility: f64,
    pub excess_max_drawdown: f64,
    pub holding_count: usize,
    pub average_weight: f64,
    pub max_weight: f64,
    pub concentration: f64,
    pub weight_std_dev: f64,
    pub median_weight: f64,
    pub average_daily_turnover: f64,
    pub total_assets: f64,
    pub cash_balance: f64,
    pub unit_nav: f64,
    pub initial_cash: f64,
    pub excess_win_rate: f64,
    pub monthly_sharpe: f64,
    pub monthly_volatility: f64,
}

pub fn compute_backtest_metrics(
    equity_curve: &[DailyEquityPoint],
    fills: &[FillEvent],
    daily_holdings: &[HoldingSummary],
    initial_cash: f64,
) -> BacktestMetrics {
    let Some(first_point) = equity_curve.first() else {
        return BacktestMetrics {
            risk_free_rate: DEFAULT_RISK_FREE_RATE,
            initial_cash,
            ..BacktestMetrics::default()
        };
    };
    let Some(last_point) = equity_curve.last() else {
        return BacktestMetrics {
            risk_free_rate: DEFAULT_RISK_FREE_RATE,
            initial_cash,
            ..BacktestMetrics::default()
        };
    };

    let trade_days = equity_curve.len();
    let benchmark_start = if first_point.benchmark_prev_close.is_finite()
        && first_point.benchmark_prev_close > f64::EPSILON
    {
        first_point.benchmark_prev_close
    } else {
        first_point.benchmark_close
    };
    let mut returns = Vec::with_capacity(equity_curve.len());
    returns.push(pct_change(initial_cash, first_point.total_equity));
    returns.extend(
        equity_curve
            .windows(2)
            .map(|window| pct_change(window[0].total_equity, window[1].total_equity)),
    );
    let mut benchmark_returns = Vec::with_capacity(equity_curve.len());
    benchmark_returns.push(pct_change(benchmark_start, first_point.benchmark_close));
    benchmark_returns.extend(
        equity_curve
            .windows(2)
            .map(|window| pct_change(window[0].benchmark_close, window[1].benchmark_close)),
    );
    let excess_returns = returns
        .iter()
        .zip(benchmark_returns.iter())
        .map(|(lhs, rhs)| lhs - rhs)
        .collect::<Vec<_>>();

    let benchmark_net_value = if benchmark_start.abs() < f64::EPSILON {
        1.0
    } else {
        last_point.benchmark_close / benchmark_start
    };
    let benchmark_cumulative_return = benchmark_net_value - 1.0;
    let total_return = if initial_cash.abs() < f64::EPSILON {
        0.0
    } else {
        (last_point.total_equity / initial_cash) - 1.0
    };
    let excess_cumulative_return = if benchmark_net_value.abs() < f64::EPSILON {
        total_return
    } else {
        (last_point.total_equity / initial_cash) / benchmark_net_value - 1.0
    };
    let excess_return = total_return - benchmark_cumulative_return;
    let annual_return = annualize_return(total_return, trade_days);
    let excess_annual_return = annualize_return(excess_cumulative_return, trade_days);

    let risk_free_rate = DEFAULT_RISK_FREE_RATE;
    let daily_rf = risk_free_rate / TRADING_DAYS_PER_YEAR;
    let sharpe = annualized_sharpe(&returns, daily_rf, TRADING_DAYS_PER_YEAR);
    let sortino = annualized_sortino(&returns, daily_rf, TRADING_DAYS_PER_YEAR);
    let information_ratio = annualized_sharpe(&excess_returns, 0.0, TRADING_DAYS_PER_YEAR);
    let tracking_error = annualized_std(&excess_returns, TRADING_DAYS_PER_YEAR);
    let volatility = annualized_std(&returns, TRADING_DAYS_PER_YEAR);
    let excess_volatility = annualized_std(&excess_returns, TRADING_DAYS_PER_YEAR);
    let excess_sharpe = annualized_sharpe(&excess_returns, 0.0, TRADING_DAYS_PER_YEAR);
    let (alpha, beta) = alpha_beta(&returns, &benchmark_returns, daily_rf);

    let equity_nav = equity_curve
        .iter()
        .map(|point| safe_div(point.total_equity, initial_cash, 1.0))
        .collect::<Vec<_>>();
    let benchmark_nav_series = equity_curve
        .iter()
        .map(|point| safe_div(point.benchmark_close, benchmark_start, 1.0))
        .collect::<Vec<_>>();
    let excess_nav_series = equity_nav
        .iter()
        .zip(benchmark_nav_series.iter())
        .map(|(lhs, rhs)| safe_div(*lhs, *rhs, *lhs))
        .collect::<Vec<_>>();

    let (max_drawdown, max_drawdown_duration_days) = drawdown_stats(&equity_nav);
    let (excess_max_drawdown, _) = drawdown_stats(&excess_nav_series);

    let winning_days = returns.iter().filter(|value| **value > 0.0).count();
    let excess_winning_days = excess_returns.iter().filter(|value| **value > 0.0).count();
    let win_rate = ratio(winning_days, returns.len());
    let excess_win_rate = ratio(excess_winning_days, excess_returns.len());

    let monthly_portfolio_returns =
        group_monthly_returns(equity_curve, initial_cash, |point| point.total_equity);
    let monthly_benchmark_returns =
        group_monthly_returns(equity_curve, benchmark_start, |point| point.benchmark_close);
    let monthly_excess_returns = monthly_portfolio_returns
        .iter()
        .zip(monthly_benchmark_returns.iter())
        .map(|(lhs, rhs)| lhs - rhs)
        .collect::<Vec<_>>();
    let monthly_excess_win_rate = ratio(
        monthly_excess_returns
            .iter()
            .filter(|value| **value > 0.0)
            .count(),
        monthly_excess_returns.len(),
    );
    let monthly_sharpe = annualized_sharpe(
        &monthly_portfolio_returns,
        risk_free_rate / MONTHS_PER_YEAR,
        MONTHS_PER_YEAR,
    );
    let monthly_volatility = annualized_std(&monthly_portfolio_returns, MONTHS_PER_YEAR);

    let turnover_by_date = fills
        .iter()
        .fold(BTreeMap::<NaiveDate, f64>::new(), |mut acc, fill| {
            *acc.entry(fill.date).or_default() += fill.gross_amount.abs();
            acc
        });
    let equity_by_date = equity_curve
        .iter()
        .map(|point| (point.date, point.total_equity))
        .collect::<BTreeMap<_, _>>();
    let average_daily_turnover = if equity_curve.is_empty() {
        0.0
    } else {
        equity_curve
            .iter()
            .map(|point| {
                let traded = turnover_by_date
                    .get(&point.date)
                    .copied()
                    .unwrap_or_default();
                safe_div(traded, point.total_equity.max(initial_cash * 0.5), 0.0)
            })
            .sum::<f64>()
            / equity_curve.len() as f64
    };

    let latest_date = last_point.date;
    let latest_holdings = daily_holdings
        .iter()
        .filter(|row| row.date == latest_date && row.quantity > 0)
        .collect::<Vec<_>>();
    let weights = latest_holdings
        .iter()
        .map(|holding| safe_div(holding.market_value, last_point.total_equity, 0.0))
        .collect::<Vec<_>>();
    let holding_count = latest_holdings.len();
    let average_weight = mean(&weights);
    let max_weight = weights
        .iter()
        .copied()
        .fold(0.0_f64, |acc, value| acc.max(value));
    let concentration = weights.iter().map(|weight| weight * weight).sum::<f64>();
    let weight_std_dev = std_dev(&weights);
    let median_weight = median(&weights);

    let total_trade_days = equity_by_date.len();

    BacktestMetrics {
        total_return,
        annual_return,
        sharpe,
        max_drawdown,
        win_rate,
        alpha,
        beta,
        benchmark_cumulative_return,
        benchmark_net_value,
        risk_free_rate,
        monthly_excess_win_rate,
        excess_cumulative_return,
        excess_annual_return,
        max_drawdown_duration_days,
        total_trade_days,
        sortino,
        information_ratio,
        tracking_error,
        volatility,
        excess_return,
        excess_sharpe,
        excess_volatility,
        excess_max_drawdown,
        holding_count,
        average_weight,
        max_weight,
        concentration,
        weight_std_dev,
        median_weight,
        average_daily_turnover,
        total_assets: last_point.total_equity,
        cash_balance: last_point.cash,
        unit_nav: safe_div(last_point.total_equity, initial_cash, 0.0),
        initial_cash,
        excess_win_rate,
        monthly_sharpe,
        monthly_volatility,
    }
}

fn pct_change(previous: f64, current: f64) -> f64 {
    if previous.abs() < f64::EPSILON {
        0.0
    } else {
        (current / previous) - 1.0
    }
}

fn annualize_return(total_return: f64, periods: usize) -> f64 {
    if periods == 0 {
        return 0.0;
    }
    let periods = periods as f64;
    let base = 1.0 + total_return;
    if base <= 0.0 {
        return -1.0;
    }
    base.powf(TRADING_DAYS_PER_YEAR / periods) - 1.0
}

fn annualized_sharpe(returns: &[f64], daily_rf: f64, periods_per_year: f64) -> f64 {
    if returns.len() < 2 {
        return 0.0;
    }
    let adjusted = returns
        .iter()
        .map(|value| value - daily_rf)
        .collect::<Vec<_>>();
    let mean_ret = mean(&adjusted);
    let std = std_dev(&adjusted);
    if std <= f64::EPSILON {
        0.0
    } else {
        mean_ret / std * periods_per_year.sqrt()
    }
}

fn annualized_sortino(returns: &[f64], daily_rf: f64, periods_per_year: f64) -> f64 {
    if returns.is_empty() {
        return 0.0;
    }
    let adjusted = returns
        .iter()
        .map(|value| value - daily_rf)
        .collect::<Vec<_>>();
    let downside = adjusted
        .iter()
        .filter(|value| **value < 0.0)
        .map(|value| value.powi(2))
        .collect::<Vec<_>>();
    if downside.is_empty() {
        return 0.0;
    }
    let downside_dev = (downside.iter().sum::<f64>() / downside.len() as f64).sqrt();
    if downside_dev <= f64::EPSILON {
        0.0
    } else {
        mean(&adjusted) / downside_dev * periods_per_year.sqrt()
    }
}

fn annualized_std(values: &[f64], periods_per_year: f64) -> f64 {
    std_dev(values) * periods_per_year.sqrt()
}

fn alpha_beta(returns: &[f64], benchmark_returns: &[f64], daily_rf: f64) -> (f64, f64) {
    if returns.len() < 2 || returns.len() != benchmark_returns.len() {
        return (0.0, 0.0);
    }
    let strategy_excess = returns
        .iter()
        .map(|value| value - daily_rf)
        .collect::<Vec<_>>();
    let benchmark_excess = benchmark_returns
        .iter()
        .map(|value| value - daily_rf)
        .collect::<Vec<_>>();
    let mean_strategy = mean(&strategy_excess);
    let mean_benchmark = mean(&benchmark_excess);
    let variance_benchmark = variance(&benchmark_excess);
    if variance_benchmark <= f64::EPSILON {
        return (0.0, 0.0);
    }
    let covariance = strategy_excess
        .iter()
        .zip(benchmark_excess.iter())
        .map(|(lhs, rhs)| (lhs - mean_strategy) * (rhs - mean_benchmark))
        .sum::<f64>()
        / (strategy_excess.len() - 1) as f64;
    let beta = covariance / variance_benchmark;
    let alpha = (mean_strategy - beta * mean_benchmark) * TRADING_DAYS_PER_YEAR;
    (alpha, beta)
}

fn drawdown_stats(nav: &[f64]) -> (f64, usize) {
    let mut peak = 0.0_f64;
    let mut max_drawdown = 0.0_f64;
    let mut duration = 0_usize;
    let mut max_duration = 0_usize;
    for value in nav {
        if *value >= peak {
            peak = *value;
            duration = 0;
            continue;
        }
        if peak > f64::EPSILON {
            let drawdown = (*value / peak) - 1.0;
            if drawdown < max_drawdown {
                max_drawdown = drawdown;
            }
        }
        duration += 1;
        if duration > max_duration {
            max_duration = duration;
        }
    }
    (max_drawdown, max_duration)
}

fn group_monthly_returns<F>(
    equity_curve: &[DailyEquityPoint],
    initial_value: f64,
    value_fn: F,
) -> Vec<f64>
where
    F: Fn(&DailyEquityPoint) -> f64,
{
    let mut month_last = BTreeMap::<(i32, u32), f64>::new();
    let mut month_first = BTreeMap::<(i32, u32), f64>::new();
    let mut previous_value = initial_value;
    for point in equity_curve {
        let key = (point.date.year(), point.date.month());
        let value = value_fn(point);
        month_first.entry(key).or_insert(previous_value);
        month_last.insert(key, value);
        previous_value = value;
    }
    let mut keys = month_last.keys().copied().collect::<Vec<_>>();
    keys.sort_unstable();
    keys.into_iter()
        .filter_map(|key| {
            let first = month_first.get(&key).copied().unwrap_or_default();
            let last = month_last.get(&key).copied().unwrap_or_default();
            if first.abs() < f64::EPSILON {
                None
            } else {
                Some((last / first) - 1.0)
            }
        })
        .collect()
}

fn mean(values: &[f64]) -> f64 {
    if values.is_empty() {
        0.0
    } else {
        values.iter().sum::<f64>() / values.len() as f64
    }
}

fn variance(values: &[f64]) -> f64 {
    if values.len() < 2 {
        return 0.0;
    }
    let avg = mean(values);
    values
        .iter()
        .map(|value| (value - avg).powi(2))
        .sum::<f64>()
        / (values.len() - 1) as f64
}

fn std_dev(values: &[f64]) -> f64 {
    variance(values).sqrt()
}

fn median(values: &[f64]) -> f64 {
    if values.is_empty() {
        return 0.0;
    }
    let mut sorted = values.to_vec();
    sorted.sort_by(|lhs, rhs| lhs.partial_cmp(rhs).unwrap_or(std::cmp::Ordering::Equal));
    let mid = sorted.len() / 2;
    if sorted.len() % 2 == 0 {
        (sorted[mid - 1] + sorted[mid]) / 2.0
    } else {
        sorted[mid]
    }
}

fn ratio(numerator: usize, denominator: usize) -> f64 {
    if denominator == 0 {
        0.0
    } else {
        numerator as f64 / denominator as f64
    }
}

fn safe_div(numerator: f64, denominator: f64, fallback: f64) -> f64 {
    if denominator.abs() < f64::EPSILON {
        fallback
    } else {
        numerator / denominator
    }
}

#[cfg(test)]
mod tests {
    use super::*;

    fn equity_point(
        date: &str,
        total_equity: f64,
        benchmark_close: f64,
        benchmark_prev_close: f64,
    ) -> DailyEquityPoint {
        DailyEquityPoint {
            date: NaiveDate::parse_from_str(date, "%Y-%m-%d").unwrap(),
            cash: total_equity,
            market_value: 0.0,
            total_equity,
            benchmark_close,
            benchmark_prev_close,
            notes: String::new(),
            diagnostics: String::new(),
        }
    }

    #[test]
    fn benchmark_cumulative_return_uses_first_day_previous_close() {
        let curve = vec![
            equity_point("2025-01-02", 100.0, 5797.089, 5957.717),
            equity_point("2025-12-31", 120.0, 7595.285, 7597.299),
        ];
        let metrics = compute_backtest_metrics(&curve, &[], &[], 100.0);
        let expected = 7595.285 / 5957.717 - 1.0;
        assert!((metrics.benchmark_cumulative_return - expected).abs() < 1e-12);
    }
}