monte carlo backtesting - why single backtest runs lie to you

from dataclasses import dataclass, field
from typing import List, Dict, Optional, Tuple
import numpy as np
import pandas as pd
from concurrent.futures import ProcessPoolExecutor
import warnings

warnings.filterwarnings('ignore')


@dataclass
class TradeResult:
    """Single trade outcome"""
    pnl: float
    duration_hours: float
    instrument: str
    entry_price: float
    exit_price: float
    position_size: float


@dataclass
class MonteCarloConfig:
    """Configuration for MC simulation"""
    n_simulations: int = 10_000
    confidence_levels: List[float] = field(
        default_factory=lambda: [0.05, 0.25, 0.50, 0.75, 0.95]
    )
    use_block_bootstrap: bool = True
    block_size: int = 5  # preserve some autocorrelation
    account_size: float = 1_200_000.0
    max_workers: int = 8


class MonteCarloBacktest:
    """Monte Carlo simulation engine for trade sequences"""

    def __init__(self, config: MonteCarloConfig = None):
        self.config = config or MonteCarloConfig()
        self.trades: List[TradeResult] = []
        self.results: Optional[Dict] = None

    def load_trades(self, trades: List[TradeResult]):
        """Load historical trade results"""
        self.trades = trades
        self._pnl_array = np.array([t.pnl for t in trades])

    def _single_simulation(self, seed: int) -> np.ndarray:
        """Run one MC path"""
        rng = np.random.RandomState(seed)
        n_trades = len(self._pnl_array)

        if self.config.use_block_bootstrap:
            # block bootstrap preserves short-term autocorrelation
            blocks = []
            while len(blocks) < n_trades:
                start = rng.randint(0, n_trades - self.config.block_size)
                block = self._pnl_array[start:start + self.config.block_size]
                blocks.extend(block)
            resampled = np.array(blocks[:n_trades])
        else:
            # simple bootstrap - iid assumption
            indices = rng.randint(0, n_trades, size=n_trades)
            resampled = self._pnl_array[indices]

        # cumulative equity curve
        equity = self.config.account_size + np.cumsum(resampled)
        return equity

    def run(self) -> Dict:
        """Run full Monte Carlo simulation"""
        seeds = range(self.config.n_simulations)

        # parallel execution
        with ProcessPoolExecutor(max_workers=self.config.max_workers) as executor:
            equity_curves = list(executor.map(self._single_simulation, seeds))

        curves = np.array(equity_curves)

        # calculate stats across all simulations
        final_values = curves[:, -1]
        returns = (final_values - self.config.account_size) / self.config.account_size

        # drawdown calculation for each path
        max_drawdowns = []
        for curve in curves:
            running_max = np.maximum.accumulate(curve)
            drawdowns = (curve - running_max) / running_max
            max_drawdowns.append(np.min(drawdowns))

        max_drawdowns = np.array(max_drawdowns)

        # percentile analysis
        percentiles = {}
        for level in self.config.confidence_levels:
            pct = int(level * 100)
            percentiles[f'p{pct}'] = {
                'final_return': float(np.percentile(returns, pct)),
                'max_drawdown': float(np.percentile(max_drawdowns, pct)),
                'final_value': float(np.percentile(final_values, pct))
            }

        # probability of ruin (losing > 25% of account)
        ruin_threshold = self.config.account_size * 0.75
        prob_ruin = np.mean(np.min(curves, axis=1) < ruin_threshold)

        # probability of profit
        prob_profit = np.mean(final_values > self.config.account_size)

        self.results = {
            'n_simulations': self.config.n_simulations,
            'n_trades': len(self.trades),
            'median_return': float(np.median(returns)),
            'mean_return': float(np.mean(returns)),
            'std_return': float(np.std(returns)),
            'median_max_drawdown': float(np.median(max_drawdowns)),
            'worst_drawdown': float(np.min(max_drawdowns)),
            'prob_profit': float(prob_profit),
            'prob_ruin': float(prob_ruin),
            'percentiles': percentiles,
            'curves': curves,  # full data for plotting
            'max_drawdowns': max_drawdowns
        }

        return self.results

    def summary(self) -> str:
        """Print human-readable summary"""
        if self.results is None:
            raise ValueError("run simulation first")

        r = self.results
        lines = [
            f"Monte Carlo Simulation ({r['n_simulations']:,} paths, {r['n_trades']} trades)",
            f"",
            f"  median return:      {r['median_return']:+.1%}",
            f"  mean return:        {r['mean_return']:+.1%}",
            f"  std of returns:     {r['std_return']:.1%}",
            f"  median max DD:      {r['median_max_drawdown']:.1%}",
            f"  worst max DD:       {r['worst_drawdown']:.1%}",
            f"  prob of profit:     {r['prob_profit']:.1%}",
            f"  prob of ruin (>25%): {r['prob_ruin']:.1%}",
            f"",
            f"  percentile breakdown:",
        ]

        for level, stats in r['percentiles'].items():
            lines.append(
                f"    {level}: return={stats['final_return']:+.1%}  "
                f"DD={stats['max_drawdown']:.1%}  "
                f"value=${stats['final_value']:,.0f}"
            )

        return "\n".join(lines)


# usage example with my actual trade log format
def load_from_csv(filepath: str) -> List[TradeResult]:
    """Load trades from my standard CSV export"""
    df = pd.read_csv(filepath)
    trades = []
    for _, row in df.iterrows():
        trades.append(TradeResult(
            pnl=row['realized_pnl'],
            duration_hours=row['hold_hours'],
            instrument=row['symbol'],
            entry_price=row['entry'],
            exit_price=row['exit'],
            position_size=row['qty']
        ))
    return trades