Merge pull request #11 from alphabench/feat/refining-metrics-for-portfolio

docs: remove iframe from README, post1 release

Cargo.lock

@@ -502,7 +502,7 @@ dependencies = [
 
 [[package]]
 name = "raptorbt"
-version = "0.3.0"
+version = "0.3.2"
 dependencies = [
  "approx",
  "criterion",

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "raptorbt"
-version = "0.3.0"
+version = "0.3.2"
 edition = "2021"
 description = "High-performance Rust backtesting engine with Python bindings. Drop-in VectorBT replacement with up insanely faster performance at fractional memory footprint."
 authors = ["Alphabench <contact@alphabench.in>"]

README.md

@@ -79,9 +79,10 @@ RaptorBT was built to address the performance limitations of VectorBT. Benchmark
 
 ### Key Features
 
-- **5 Strategy Types**: Single instrument, basket/collective, pairs trading, options, and multi-strategy
-- **30+ Metrics**: Full parity with VectorBT including Sharpe, Sortino, Calmar, Omega, SQN, and more
-- **10 Technical Indicators**: SMA, EMA, RSI, MACD, Stochastic, ATR, Bollinger Bands, ADX, VWAP, Supertrend
+- **6 Strategy Types**: Single instrument, basket/collective, pairs trading, options, spreads, and multi-strategy
+- **Monte Carlo Simulation**: Correlated multi-asset forward projection via GBM + Cholesky decomposition
+- **33 Metrics**: Full parity with VectorBT including Sharpe, Sortino, Calmar, Omega, SQN, Payoff Ratio, Recovery Factor, and more
+- **12 Technical Indicators**: SMA, EMA, RSI, MACD, Stochastic, ATR, Bollinger Bands, ADX, VWAP, Supertrend, Rolling Min, Rolling Max
 - **Stop/Target Management**: Fixed, ATR-based, and trailing stops with risk-reward targets
 - **100% Deterministic**: No JIT compilation variance between runs
 - **Native Parallelism**: Rayon-based parallel processing with explicit SIMD optimizations
@@ -127,6 +128,7 @@ raptorbt/
 │   ├── core/              # Core types and error handling
 │   │   ├── types.rs       # BacktestConfig, BacktestResult, Trade, Metrics
 │   │   ├── error.rs       # RaptorError enum
+│   │   ├── session.rs     # SessionTracker, SessionConfig (intraday sessions)
 │   │   └── timeseries.rs  # Time series utilities
 │   │
 │   ├── strategies/        # Strategy implementations
@@ -134,6 +136,7 @@ raptorbt/
 │   │   ├── basket.rs      # Basket/collective strategies
 │   │   ├── pairs.rs       # Pairs trading
 │   │   ├── options.rs     # Options strategies
+│   │   ├── spreads.rs     # Multi-leg spread strategies
 │   │   └── multi.rs       # Multi-strategy combining
 │   │
 │   ├── indicators/        # Technical indicators
@@ -141,7 +144,8 @@ raptorbt/
 │   │   ├── momentum.rs    # RSI, MACD, Stochastic
 │   │   ├── volatility.rs  # ATR, Bollinger Bands
 │   │   ├── strength.rs    # ADX
-│   │   └── volume.rs      # VWAP
+│   │   ├── volume.rs      # VWAP
+│   │   └── rolling.rs     # Rolling Min/Max (LLV/HHV)
 │   │
 │   ├── metrics/           # Performance metrics
 │   │   ├── streaming.rs   # Streaming metric calculations
@@ -158,6 +162,12 @@ raptorbt/
 │   │   ├── atr.rs         # ATR-based stops
 │   │   └── trailing.rs    # Trailing stops
 │   │
+│   ├── portfolio/         # Portfolio-level analysis
+│   │   ├── monte_carlo.rs # Monte Carlo forward simulation (GBM + Cholesky)
+│   │   ├── allocation.rs  # Capital allocation
+│   │   ├── engine.rs      # Portfolio engine
+│   │   └── position.rs    # Position management
+│   │
 │   ├── python/            # PyO3 bindings
 │   │   ├── bindings.rs    # Python function exports
 │   │   └── numpy_bridge.rs # NumPy array conversion
@@ -529,6 +539,68 @@ config.set_risk_reward_target(ratio=2.0)  # 2:1 risk-reward ratio
 
 ---
 
+## Monte Carlo Portfolio Simulation
+
+RaptorBT includes a high-performance Monte Carlo forward simulation engine for portfolio risk analysis. It uses Geometric Brownian Motion (GBM) with Cholesky decomposition for correlated multi-asset simulation, parallelized via Rayon.
+
+```python
+import numpy as np
+import raptorbt
+
+# Historical daily returns per strategy/asset (numpy arrays)
+returns = [
+    np.array([0.001, -0.002, 0.003, ...]),  # Strategy 1 returns
+    np.array([0.002, 0.001, -0.001, ...]),   # Strategy 2 returns
+]
+
+# Portfolio weights (must sum to 1.0)
+weights = np.array([0.6, 0.4])
+
+# Correlation matrix (N x N)
+correlation_matrix = [
+    np.array([1.0, 0.3]),
+    np.array([0.3, 1.0]),
+]
+
+# Run simulation
+result = raptorbt.simulate_portfolio_mc(
+    returns=returns,
+    weights=weights,
+    correlation_matrix=correlation_matrix,
+    initial_value=100000.0,
+    n_simulations=10000,   # Number of Monte Carlo paths (default: 10,000)
+    horizon_days=252,      # Forward projection horizon (default: 252)
+    seed=42,               # Random seed for reproducibility (default: 42)
+)
+
+# Results
+print(f"Expected Return: {result['expected_return']:.2f}%")
+print(f"Probability of Loss: {result['probability_of_loss']:.2%}")
+print(f"VaR (95%): {result['var_95']:.2f}%")
+print(f"CVaR (95%): {result['cvar_95']:.2f}%")
+
+# Percentile paths: list of (percentile, path_values)
+# Percentiles: 5th, 25th, 50th, 75th, 95th
+for pct, path in result['percentile_paths']:
+    print(f"  P{pct:.0f} final value: {path[-1]:.2f}")
+
+# Final values: numpy array of terminal values for all simulations
+final_values = result['final_values']  # numpy array, length = n_simulations
+```
+
+### Result Fields
+
+| Field                 | Type                       | Description                                                |
+| --------------------- | -------------------------- | ---------------------------------------------------------- |
+| `expected_return`     | `float`                    | Expected return as percentage over the horizon             |
+| `probability_of_loss` | `float`                    | Probability that final value < initial value (0.0 to 1.0)  |
+| `var_95`              | `float`                    | Value at Risk at 95% confidence (percentage)               |
+| `cvar_95`             | `float`                    | Conditional VaR at 95% confidence (percentage)             |
+| `percentile_paths`    | `List[Tuple[float, List]]` | Portfolio paths at 5th, 25th, 50th, 75th, 95th percentiles |
+| `final_values`        | `numpy.ndarray`            | Terminal portfolio values for all simulations              |
+
+---
+
 ## VectorBT Comparison
 
 RaptorBT is designed as a drop-in replacement for VectorBT. Here's a side-by-side comparison:
@@ -643,10 +715,27 @@ inst_config.set_fixed_target(0.05)
 ```
 
 **Fields:**
+
 - `lot_size` - Minimum tradeable quantity. Position sizes are rounded down to nearest lot_size multiple. Use `1.0` for equities, `50.0` for NIFTY F&O, `0.01` for forex.
 - `alloted_capital` - Per-instrument capital cap (capped at available cash).
 - `existing_qty` / `avg_price` - Reserved for future live-to-backtest transitions.
 
+### simulate_portfolio_mc
+
+```python
+result = raptorbt.simulate_portfolio_mc(
+    returns: List[np.ndarray],               # Per-asset daily returns (N arrays)
+    weights: np.ndarray,                     # Portfolio weights (length N, sum to 1)
+    correlation_matrix: List[np.ndarray],    # N x N correlation matrix
+    initial_value: float,                    # Starting portfolio value
+    n_simulations: int = 10000,              # Number of Monte Carlo paths
+    horizon_days: int = 252,                 # Forward projection horizon in days
+    seed: int = 42,                          # Random seed for reproducibility
+) -> dict
+```
+
+Returns a dictionary with keys: `expected_return`, `probability_of_loss`, `var_95`, `cvar_95`, `percentile_paths`, `final_values`.
+
 ### PyBacktestResult
 
 ```python
@@ -699,6 +788,8 @@ metrics.max_consecutive_wins
 metrics.max_consecutive_losses
 metrics.exposure_pct
 metrics.open_trade_pnl
+metrics.payoff_ratio            # avg win / avg loss (risk/reward per trade)
+metrics.recovery_factor         # net profit / max drawdown (resilience)
 
 # Convert to dictionary (VectorBT format)
 stats_dict = metrics.to_dict()
@@ -833,6 +924,22 @@ MIT License - see [LICENSE](LICENSE) for details.
 
 ## Changelog
 
+### v0.3.2
+
+- Add `payoff_ratio` metric to `BacktestMetrics` — average winning trade return divided by average losing trade return (absolute), measures risk/reward per trade
+- Add `recovery_factor` metric to `BacktestMetrics` — net profit divided by maximum drawdown in absolute terms, measures how many times over the strategy recovered from its worst drawdown
+- Both metrics computed in `StreamingMetrics::finalize()` (single-instrument backtest) and `PortfolioEngine` (multi-strategy aggregation)
+- Both metrics exposed via PyO3 as `#[pyo3(get)]` attributes on `PyBacktestMetrics`
+- Handles edge cases: returns `f64::INFINITY` when denominator is zero with positive numerator, `0.0` otherwise
+
+### v0.3.1
+
+- Add Monte Carlo portfolio simulation (`simulate_portfolio_mc`) for forward risk projection
+- Geometric Brownian Motion (GBM) with Cholesky decomposition for correlated multi-asset simulation
+- Rayon-parallelized simulation paths with deterministic seeding (xoshiro256\*\*)
+- Returns percentile paths (P5/P25/P50/P75/P95), VaR, CVaR, expected return, and probability of loss
+- GIL released during simulation for maximum Python concurrency
+
 ### v0.3.0
 
 - Per-instrument configuration via `PyInstrumentConfig` (lot_size, alloted_capital, stop/target overrides)

pyproject.toml

@@ -4,7 +4,7 @@ build-backend = "maturin"
 
 [project]
 name = "raptorbt"
-version = "0.3.0"
+version = "0.3.2.post1"
 description = "High-performance Rust backtesting engine with Python bindings. Drop-in VectorBT replacement with up insanely faster performance at fractional memory footprint."
 readme = "README.md"
 requires-python = ">=3.10"

python/raptorbt/__init__.py

@@ -26,6 +26,8 @@ from raptorbt._raptorbt import (
     run_pairs_backtest,
     run_multi_backtest,
     run_spread_backtest,
+    # Monte Carlo simulation
+    simulate_portfolio_mc,
     # Indicator functions
     sma,
     ema,
@@ -41,7 +43,7 @@ from raptorbt._raptorbt import (
     rolling_max,
 )
 
-__version__ = "0.3.0"
+__version__ = "0.3.2.post1"
 
 __all__ = [
     # Config classes
@@ -60,6 +62,8 @@ __all__ = [
     "run_pairs_backtest",
     "run_multi_backtest",
     "run_spread_backtest",
+    # Monte Carlo simulation
+    "simulate_portfolio_mc",
     # Indicator functions
     "sma",
     "ema",

src/core/types.rs

@@ -376,6 +376,10 @@ pub struct BacktestMetrics {
     pub avg_holding_period: f64,
     /// Exposure time percentage (time in market).
     pub exposure_pct: f64,
+    /// Payoff ratio (avg win / avg loss).
+    pub payoff_ratio: f64,
+    /// Recovery factor (net profit / max drawdown).
+    pub recovery_factor: f64,
 }
 
 /// Complete backtest result.

src/lib.rs

@@ -44,6 +44,9 @@ fn _raptorbt(_py: Python<'_>, m: &PyModule) -> PyResult<()> {
     m.add_function(wrap_pyfunction!(python::bindings::run_multi_backtest, m)?)?;
     m.add_function(wrap_pyfunction!(python::bindings::run_spread_backtest, m)?)?;
 
+    // Register Monte Carlo simulation
+    m.add_function(wrap_pyfunction!(python::bindings::simulate_portfolio_mc, m)?)?;
+
     // Register indicator functions
     m.add_function(wrap_pyfunction!(python::bindings::sma, m)?)?;
     m.add_function(wrap_pyfunction!(python::bindings::ema, m)?)?;

src/metrics/streaming.rs

@@ -513,6 +513,30 @@ impl StreamingMetrics {
         let worst_trade_pct =
             if self.worst_trade_pct == f64::INFINITY { 0.0 } else { self.worst_trade_pct };
 
+        // Payoff ratio: average win / average loss (absolute value)
+        let payoff_ratio = if avg_loss_pct.abs() > 0.0 {
+            avg_win_pct / avg_loss_pct.abs()
+        } else if avg_win_pct > 0.0 {
+            f64::INFINITY
+        } else {
+            0.0
+        };
+
+        // Recovery factor: net profit / max drawdown (absolute value)
+        let net_profit = final_value - initial_capital;
+        let recovery_factor = if self.max_drawdown_pct > 0.0 && initial_capital > 0.0 {
+            let max_dd_absolute = self.max_drawdown_pct / 100.0 * initial_capital;
+            if max_dd_absolute > 0.0 {
+                net_profit / max_dd_absolute
+            } else {
+                0.0
+            }
+        } else if net_profit > 0.0 {
+            f64::INFINITY
+        } else {
+            0.0
+        };
+
         BacktestMetrics {
             total_return_pct,
             sharpe_ratio,
@@ -545,6 +569,8 @@ impl StreamingMetrics {
             max_consecutive_losses: self.max_consecutive_losses,
             avg_holding_period,
             exposure_pct: 0.0, // TODO: calculate based on time in market
+            payoff_ratio,
+            recovery_factor,
         }
     }
 

src/portfolio/engine.rs

@@ -591,6 +591,30 @@ impl PortfolioEngine {
             0.0
         };
 
+        // Payoff ratio: average win / average loss (absolute value)
+        let payoff_ratio = if avg_loss_pct.abs() > 0.0 {
+            avg_win_pct / avg_loss_pct.abs()
+        } else if avg_win_pct > 0.0 {
+            f64::INFINITY
+        } else {
+            0.0
+        };
+
+        // Recovery factor: net profit / max drawdown (absolute value)
+        let net_profit = end_value - start_value;
+        let recovery_factor = if max_drawdown_pct > 0.0 && start_value > 0.0 {
+            let max_dd_absolute = max_drawdown_pct / 100.0 * start_value;
+            if max_dd_absolute > 0.0 {
+                net_profit / max_dd_absolute
+            } else {
+                0.0
+            }
+        } else if net_profit > 0.0 {
+            f64::INFINITY
+        } else {
+            0.0
+        };
+
         BacktestMetrics {
             total_return_pct,
             sharpe_ratio,
@@ -623,6 +647,8 @@ impl PortfolioEngine {
             max_consecutive_losses,
             avg_holding_period,
             exposure_pct,
+            payoff_ratio,
+            recovery_factor,
         }
     }
 

src/portfolio/mod.rs

@@ -2,8 +2,10 @@
 
 pub mod allocation;
 pub mod engine;
+pub mod monte_carlo;
 pub mod position;
 
 pub use allocation::{AllocationStrategy, CapitalAllocator};
 pub use engine::PortfolioEngine;
+pub use monte_carlo::{simulate_portfolio_forward, MonteCarloConfig, MonteCarloResult};
 pub use position::PositionManager;

src/portfolio/monte_carlo.rs

@@ -0,0 +1,361 @@
+//! Monte Carlo forward simulation for portfolio projection.
+//!
+//! Uses Geometric Brownian Motion (GBM) with Cholesky decomposition
+//! for correlated multi-asset simulation. Parallelized via Rayon.
+
+use rayon::prelude::*;
+
+/// Configuration for Monte Carlo simulation.
+#[derive(Debug, Clone)]
+pub struct MonteCarloConfig {
+    pub n_simulations: usize,
+    pub horizon_days: usize,
+    pub seed: u64,
+}
+
+impl Default for MonteCarloConfig {
+    fn default() -> Self {
+        Self { n_simulations: 10_000, horizon_days: 252, seed: 42 }
+    }
+}
+
+/// Result of a Monte Carlo simulation.
+#[derive(Debug, Clone)]
+pub struct MonteCarloResult {
+    /// Percentile paths: Vec of (percentile, path_values)
+    pub percentile_paths: Vec<(f64, Vec<f64>)>,
+    /// Terminal value for each simulation
+    pub final_values: Vec<f64>,
+    /// Expected annualized return
+    pub expected_return: f64,
+    /// Probability of loss (final value < initial value)
+    pub probability_of_loss: f64,
+    /// Value at Risk at 95% confidence
+    pub var_95: f64,
+    /// Conditional Value at Risk at 95% confidence
+    pub cvar_95: f64,
+}
+
+/// Cholesky decomposition of a symmetric positive-definite matrix.
+/// Returns lower-triangular matrix L such that A = L * L^T.
+fn cholesky(matrix: &[Vec<f64>]) -> Result<Vec<Vec<f64>>, &'static str> {
+    let n = matrix.len();
+    let mut l = vec![vec![0.0; n]; n];
+
+    for i in 0..n {
+        for j in 0..=i {
+            let mut sum = 0.0;
+            for k in 0..j {
+                sum += l[i][k] * l[j][k];
+            }
+
+            if i == j {
+                let diag = matrix[i][i] - sum;
+                if diag <= 0.0 {
+                    // Matrix is not positive definite; use a small epsilon
+                    l[i][j] = (diag.abs().max(1e-10)).sqrt();
+                } else {
+                    l[i][j] = diag.sqrt();
+                }
+            } else {
+                if l[j][j].abs() < 1e-15 {
+                    l[i][j] = 0.0;
+                } else {
+                    l[i][j] = (matrix[i][j] - sum) / l[j][j];
+                }
+            }
+        }
+    }
+
+    Ok(l)
+}
+
+/// Simple xoshiro256** PRNG for deterministic parallel simulation.
+#[derive(Clone)]
+struct Xoshiro256 {
+    s: [u64; 4],
+}
+
+impl Xoshiro256 {
+    fn new(seed: u64) -> Self {
+        // SplitMix64 to seed all 4 state words
+        let mut z = seed;
+        let mut s = [0u64; 4];
+        for item in &mut s {
+            z = z.wrapping_add(0x9e3779b97f4a7c15);
+            z = (z ^ (z >> 30)).wrapping_mul(0xbf58476d1ce4e5b9);
+            z = (z ^ (z >> 27)).wrapping_mul(0x94d049bb133111eb);
+            *item = z ^ (z >> 31);
+        }
+        Self { s }
+    }
+
+    fn jump(&mut self) {
+        // Jump function: advances state by 2^128 calls
+        const JUMP: [u64; 4] =
+            [0x180ec6d33cfd0aba, 0xd5a61266f0c9392c, 0xa9582618e03fc9aa, 0x39abdc4529b1661c];
+        let mut s0: u64 = 0;
+        let mut s1: u64 = 0;
+        let mut s2: u64 = 0;
+        let mut s3: u64 = 0;
+        for j in &JUMP {
+            for b in 0..64 {
+                if j & (1u64 << b) != 0 {
+                    s0 ^= self.s[0];
+                    s1 ^= self.s[1];
+                    s2 ^= self.s[2];
+                    s3 ^= self.s[3];
+                }
+                self.next_u64();
+            }
+        }
+        self.s[0] = s0;
+        self.s[1] = s1;
+        self.s[2] = s2;
+        self.s[3] = s3;
+    }
+
+    fn next_u64(&mut self) -> u64 {
+        let result = (self.s[1].wrapping_mul(5)).rotate_left(7).wrapping_mul(9);
+        let t = self.s[1] << 17;
+        self.s[2] ^= self.s[0];
+        self.s[3] ^= self.s[1];
+        self.s[1] ^= self.s[2];
+        self.s[0] ^= self.s[3];
+        self.s[2] ^= t;
+        self.s[3] = self.s[3].rotate_left(45);
+        result
+    }
+
+    /// Generate uniform f64 in [0, 1).
+    fn next_f64(&mut self) -> f64 {
+        (self.next_u64() >> 11) as f64 * (1.0 / (1u64 << 53) as f64)
+    }
+
+    /// Box-Muller transform for standard normal.
+    fn next_normal(&mut self) -> f64 {
+        let u1 = self.next_f64().max(1e-15);
+        let u2 = self.next_f64();
+        (-2.0 * u1.ln()).sqrt() * (2.0 * std::f64::consts::PI * u2).cos()
+    }
+}
+
+/// Core Monte Carlo simulation function.
+///
+/// # Arguments
+/// * `returns` - Per-strategy daily returns (N strategies x T days each)
+/// * `weights` - Portfolio weights (length N, must sum to 1)
+/// * `correlation_matrix` - N x N correlation matrix
+/// * `initial_value` - Starting portfolio value
+/// * `config` - Simulation configuration
+pub fn simulate_portfolio_forward(
+    returns: &[Vec<f64>],
+    weights: &[f64],
+    correlation_matrix: &[Vec<f64>],
+    initial_value: f64,
+    config: &MonteCarloConfig,
+) -> MonteCarloResult {
+    let n_assets = returns.len();
+    let dt = 1.0; // daily time step
+
+    // Compute per-asset mean and std of historical returns
+    let mut mus = vec![0.0; n_assets];
+    let mut sigmas = vec![0.0; n_assets];
+    for (i, ret) in returns.iter().enumerate() {
+        if ret.is_empty() {
+            continue;
+        }
+        let mean = ret.iter().sum::<f64>() / ret.len() as f64;
+        let var = ret.iter().map(|r| (r - mean).powi(2)).sum::<f64>() / ret.len() as f64;
+        mus[i] = mean;
+        sigmas[i] = var.sqrt().max(1e-10);
+    }
+
+    // Cholesky decomposition of correlation matrix
+    let chol = cholesky(correlation_matrix).unwrap_or_else(|_| {
+        // Fallback: identity matrix (independent assets)
+        let mut identity = vec![vec![0.0; n_assets]; n_assets];
+        for i in 0..n_assets {
+            identity[i][i] = 1.0;
+        }
+        identity
+    });
+
+    // Prepare a base RNG and create per-chunk seeds via jumping
+    let mut base_rng = Xoshiro256::new(config.seed);
+    let n_chunks = rayon::current_num_threads().max(1);
+    let chunk_size = (config.n_simulations + n_chunks - 1) / n_chunks;
+
+    let chunk_rngs: Vec<Xoshiro256> = (0..n_chunks)
+        .map(|_| {
+            let rng = base_rng.clone();
+            base_rng.jump();
+            rng
+        })
+        .collect();
+
+    // Run simulations in parallel chunks
+    let all_paths: Vec<Vec<f64>> = chunk_rngs
+        .into_par_iter()
+        .enumerate()
+        .flat_map(|(chunk_idx, mut rng)| {
+            let start = chunk_idx * chunk_size;
+            let end = (start + chunk_size).min(config.n_simulations);
+            let mut chunk_paths = Vec::with_capacity(end - start);
+
+            for _ in start..end {
+                let mut portfolio_value = initial_value;
+                let mut path = Vec::with_capacity(config.horizon_days + 1);
+                path.push(portfolio_value);
+
+                for _ in 0..config.horizon_days {
+                    // Generate N independent standard normals
+                    let z_indep: Vec<f64> = (0..n_assets).map(|_| rng.next_normal()).collect();
+
+                    // Correlate via Cholesky: z_corr = L * z_indep
+                    let mut z_corr = vec![0.0; n_assets];
+                    for i in 0..n_assets {
+                        for j in 0..=i {
+                            z_corr[i] += chol[i][j] * z_indep[j];
+                        }
+                    }
+
+                    // GBM per asset, then weighted portfolio return
+                    let mut portfolio_return = 0.0;
+                    for i in 0..n_assets {
+                        let drift = (mus[i] - 0.5 * sigmas[i].powi(2)) * dt;
+                        let diffusion = sigmas[i] * dt.sqrt() * z_corr[i];
+                        let asset_return = (drift + diffusion).exp() - 1.0;
+                        portfolio_return += weights[i] * asset_return;
+                    }
+
+                    portfolio_value *= 1.0 + portfolio_return;
+                    path.push(portfolio_value);
+                }
+
+                chunk_paths.push(path);
+            }
+
+            chunk_paths
+        })
+        .collect();
+
+    // Extract final values
+    let mut final_values: Vec<f64> = all_paths.iter().map(|p| *p.last().unwrap()).collect();
+    final_values.sort_by(|a, b| a.partial_cmp(b).unwrap_or(std::cmp::Ordering::Equal));
+
+    let n = final_values.len();
+
+    // Percentile paths: find simulations closest to each percentile's final value
+    let percentiles = [5.0, 25.0, 50.0, 75.0, 95.0];
+    let percentile_paths: Vec<(f64, Vec<f64>)> = percentiles
+        .iter()
+        .map(|&pct| {
+            let idx = ((pct / 100.0) * (n as f64 - 1.0)).round() as usize;
+            let target_final = final_values[idx.min(n - 1)];
+
+            // Find the simulation path whose final value is closest to target
+            let best_idx = all_paths
+                .iter()
+                .enumerate()
+                .min_by(|(_, a), (_, b)| {
+                    let da = (a.last().unwrap() - target_final).abs();
+                    let db = (b.last().unwrap() - target_final).abs();
+                    da.partial_cmp(&db).unwrap_or(std::cmp::Ordering::Equal)
+                })
+                .map(|(i, _)| i)
+                .unwrap_or(0);
+
+            (pct, all_paths[best_idx].clone())
+        })
+        .collect();
+
+    // Expected return (annualized from mean of final values)
+    let mean_final = final_values.iter().sum::<f64>() / n as f64;
+    let expected_return = (mean_final / initial_value - 1.0) * 100.0;
+
+    // Probability of loss
+    let n_loss = final_values.iter().filter(|&&v| v < initial_value).count();
+    let probability_of_loss = n_loss as f64 / n as f64;
+
+    // VaR 95%: 5th percentile loss
+    let p5_idx = ((0.05 * (n as f64 - 1.0)).round() as usize).min(n - 1);
+    let var_95 = ((initial_value - final_values[p5_idx]) / initial_value * 100.0).max(0.0);
+
+    // CVaR 95%: average of losses below VaR
+    let cvar_values = &final_values[..=p5_idx];
+    let cvar_95 = if cvar_values.is_empty() {
+        var_95
+    } else {
+        let avg_tail = cvar_values.iter().sum::<f64>() / cvar_values.len() as f64;
+        ((initial_value - avg_tail) / initial_value * 100.0).max(0.0)
+    };
+
+    MonteCarloResult {
+        percentile_paths,
+        final_values,
+        expected_return,
+        probability_of_loss,
+        var_95,
+        cvar_95,
+    }
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    #[test]
+    fn test_cholesky_identity() {
+        let matrix = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
+        let l = cholesky(&matrix).unwrap();
+        assert!((l[0][0] - 1.0).abs() < 1e-10);
+        assert!((l[1][1] - 1.0).abs() < 1e-10);
+        assert!(l[0][1].abs() < 1e-10);
+        assert!(l[1][0].abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_cholesky_correlated() {
+        let matrix = vec![vec![1.0, 0.5], vec![0.5, 1.0]];
+        let l = cholesky(&matrix).unwrap();
+        // Verify L * L^T = matrix
+        let reconstructed_00 = l[0][0] * l[0][0];
+        let reconstructed_01 = l[1][0] * l[0][0];
+        let reconstructed_11 = l[1][0] * l[1][0] + l[1][1] * l[1][1];
+        assert!((reconstructed_00 - 1.0).abs() < 1e-10);
+        assert!((reconstructed_01 - 0.5).abs() < 1e-10);
+        assert!((reconstructed_11 - 1.0).abs() < 1e-10);
+    }
+
+    #[test]
+    fn test_simulate_basic() {
+        // Two assets with identical positive returns
+        let returns = vec![vec![0.001; 252], vec![0.001; 252]];
+        let weights = vec![0.5, 0.5];
+        let corr = vec![vec![1.0, 0.0], vec![0.0, 1.0]];
+        let config = MonteCarloConfig { n_simulations: 100, horizon_days: 10, seed: 42 };
+
+        let result = simulate_portfolio_forward(&returns, &weights, &corr, 100000.0, &config);
+
+        assert_eq!(result.final_values.len(), 100);
+        assert_eq!(result.percentile_paths.len(), 5);
+        // Expected return should be positive given positive drift
+        assert!(result.expected_return > -50.0); // Sanity check
+    }
+
+    #[test]
+    fn test_deterministic() {
+        let returns = vec![vec![0.001; 100], vec![-0.0005; 100]];
+        let weights = vec![0.6, 0.4];
+        let corr = vec![vec![1.0, -0.3], vec![-0.3, 1.0]];
+        let config = MonteCarloConfig { n_simulations: 50, horizon_days: 20, seed: 123 };
+
+        let r1 = simulate_portfolio_forward(&returns, &weights, &corr, 100000.0, &config);
+        let r2 = simulate_portfolio_forward(&returns, &weights, &corr, 100000.0, &config);
+
+        // Same seed should produce same final values (single-threaded determinism)
+        // Note: with rayon, parallelism may affect order but not values
+        assert!((r1.expected_return - r2.expected_return).abs() < 1e-6);
+    }
+}

src/python/bindings.rs

@@ -419,6 +419,10 @@ pub struct PyBacktestMetrics {
     pub avg_holding_period: f64,
     #[pyo3(get)]
     pub exposure_pct: f64,
+    #[pyo3(get)]
+    pub payoff_ratio: f64,
+    #[pyo3(get)]
+    pub recovery_factor: f64,
 }
 
 #[pymethods]
@@ -1098,6 +1102,77 @@ pub fn rolling_max<'py>(
 // Helper Functions
 // ============================================================================
 
+// ============================================================================
+// Monte Carlo Forward Simulation
+// ============================================================================
+
+/// Run Monte Carlo forward simulation for a portfolio.
+///
+/// Uses Geometric Brownian Motion with Cholesky-decomposed correlated random
+/// draws, parallelized via Rayon.
+///
+/// # Arguments
+/// * `returns` - List of per-strategy return arrays (N strategies)
+/// * `weights` - Portfolio weight vector (length N, sums to 1)
+/// * `correlation_matrix` - N x N correlation matrix (flattened row-major as 2D list)
+/// * `initial_value` - Starting portfolio value
+/// * `n_simulations` - Number of simulation paths (default: 10000)
+/// * `horizon_days` - Forward simulation horizon in trading days (default: 252)
+/// * `seed` - Random seed for reproducibility (default: 42)
+#[pyfunction]
+#[pyo3(signature = (returns, weights, correlation_matrix, initial_value, n_simulations=10000, horizon_days=252, seed=42))]
+pub fn simulate_portfolio_mc(
+    py: Python<'_>,
+    returns: Vec<PyReadonlyArray1<'_, f64>>,
+    weights: PyReadonlyArray1<'_, f64>,
+    correlation_matrix: Vec<PyReadonlyArray1<'_, f64>>,
+    initial_value: f64,
+    n_simulations: usize,
+    horizon_days: usize,
+    seed: u64,
+) -> PyResult<PyObject> {
+    use crate::portfolio::monte_carlo::{simulate_portfolio_forward, MonteCarloConfig};
+
+    // Convert numpy arrays to Rust vecs
+    let rust_returns: Vec<Vec<f64>> =
+        returns.iter().map(|arr| arr.as_slice().unwrap().to_vec()).collect();
+
+    let rust_weights: Vec<f64> = weights.as_slice().unwrap().to_vec();
+
+    let rust_corr: Vec<Vec<f64>> =
+        correlation_matrix.iter().map(|arr| arr.as_slice().unwrap().to_vec()).collect();
+
+    let config = MonteCarloConfig { n_simulations, horizon_days, seed };
+
+    // Run simulation (releases GIL for Rayon parallelism)
+    let result = py.allow_threads(|| {
+        simulate_portfolio_forward(&rust_returns, &rust_weights, &rust_corr, initial_value, &config)
+    });
+
+    // Build Python dict result
+    let dict = pyo3::types::PyDict::new(py);
+
+    // percentile_paths: list of (percentile, list[float])
+    let paths_list = pyo3::types::PyList::empty(py);
+    for (pct, path) in &result.percentile_paths {
+        let path_list = pyo3::types::PyList::new(py, path);
+        let tuple = pyo3::types::PyTuple::new(py, &[pct.to_object(py), path_list.to_object(py)]);
+        paths_list.append(tuple)?;
+    }
+    dict.set_item("percentile_paths", paths_list)?;
+
+    // final_values as numpy array for efficiency
+    let final_arr = PyArray1::from_vec(py, result.final_values);
+    dict.set_item("final_values", final_arr)?;
+
+    dict.set_item("expected_return", result.expected_return)?;
+    dict.set_item("probability_of_loss", result.probability_of_loss)?;
+    dict.set_item("var_95", result.var_95)?;
+    dict.set_item("cvar_95", result.cvar_95)?;
+
+    Ok(dict.into())
+}
+
 /// Convert Rust BacktestResult to Python PyBacktestResult.
 fn convert_result(result: crate::core::types::BacktestResult) -> PyBacktestResult {
     let metrics = PyBacktestMetrics {
@@ -1132,6 +1207,8 @@ fn convert_result(result: crate::core::types::BacktestResult) -> PyBacktestResul
         max_consecutive_losses: result.metrics.max_consecutive_losses,
         avg_holding_period: result.metrics.avg_holding_period,
         exposure_pct: result.metrics.exposure_pct,
+        payoff_ratio: result.metrics.payoff_ratio,
+        recovery_factor: result.metrics.recovery_factor,
     };
 
     let trades: Vec<PyTrade> = result