release: v0.3.0

examples/01_trend_following.py

@@ -42,13 +42,14 @@ config = mbt.BacktestConfig(
     universe=[1],
     time_range_start=start,
     time_range_end=end,
-    bar_interval=Interval.hours(12),
+    bar_interval=Interval.hours(1),
     initial_capital=10_000,
     execution=mbt.ExecutionConfig(
         allow_short=False,
         max_position_pct=0.5,
         position_sizing_mode="FractionOfInitialCapital",
     ),
+    output_resolution=Interval.hours(1),
     fees=mbt.FeeConfig.binance_perps(),
     slippage=Slippage.fixed_bps(2),
     warmup_bars=30,

examples/13_stochastic_simulation.py

@@ -0,0 +1,145 @@
+"""Stochastic Simulation -- synthetic price paths via SDE expression DSL.
+
+Demonstrates:
+  - Built-in presets (GBM, Heston, Merton, GARCH-JD)
+  - Custom SDE model via string expressions
+  - Stochastic fan chart visualization
+  - CUDA GPU acceleration (device="cuda")
+  - All expressions compile to native Rust — full Rayon / CUDA parallelism
+
+Usage:
+    python examples/13_stochastic_simulation.py
+"""
+import time
+import manifoldbt as mbt
+
+N = 10_000_000  # 10M paths
+DEVICE = "cuda"  # "cpu" or "cuda"
+
+if __name__ == "__main__":
+    # ── 1. Geometric Brownian Motion (preset) ───────────────────────────────
+    print(f"1. GBM preset ({N:,} paths, 252 steps) [{DEVICE}]")
+    t0 = time.perf_counter()
+    result = mbt.run_stochastic(
+        "gbm",
+        s0=100.0,
+        n_paths=N,
+        n_steps=252,
+        dt=1 / 252,
+        params={"mu": 0.05, "sigma": 0.20},
+        seed=42,
+        device=DEVICE,
+    )
+    elapsed = time.perf_counter() - t0
+    print(f"   Mean final price: {result['final_price']['mean']:.2f}")
+    print(f"   Median max DD:    {result['max_drawdown']['percentiles'][3][1]:.2%}")
+    print(f"   Elapsed: {elapsed:.3f}s\n")
+
+    # ── 2. Heston stochastic volatility (preset) ────────────────────────────
+    print(f"2. Heston preset ({N:,} paths) [{DEVICE}]")
+    t0 = time.perf_counter()
+    result = mbt.run_stochastic(
+        "heston",
+        s0=100.0,
+        n_paths=N,
+        n_steps=252,
+        dt=1 / 252,
+        params={"mu": 0.05, "kappa": 2.0, "theta": 0.04, "xi": 0.3},
+        seed=42,
+        device=DEVICE,
+    )
+    elapsed = time.perf_counter() - t0
+    print(f"   Mean final price: {result['final_price']['mean']:.2f}")
+    print(f"   Ann. vol (mean):  {result['annualized_vol']['mean']:.2%}")
+    print(f"   Elapsed: {elapsed:.3f}s\n")
+
+    # ── 3. Merton Jump Diffusion (preset) ───────────────────────────────────
+    print(f"3. Merton Jump Diffusion ({N:,} paths) [{DEVICE}]")
+    t0 = time.perf_counter()
+    result = mbt.run_stochastic(
+        "merton",
+        s0=100.0,
+        n_paths=N,
+        n_steps=252,
+        dt=1 / 252,
+        params={
+            "mu": 0.05,
+            "sigma": 0.20,
+            "lambda": 1.0,     # 1 jump/year on average
+            "mu_j": -0.05,     # mean jump = -5%
+            "sigma_j": 0.08,   # jump vol = 8%
+        },
+        seed=42,
+        device=DEVICE,
+    )
+    elapsed = time.perf_counter() - t0
+    print(f"   Mean final price: {result['final_price']['mean']:.2f}")
+    print(f"   Elapsed: {elapsed:.3f}s\n")
+
+    # ── 4. Custom GARCH(1,1) Jump Diffusion ─────────────────────────────────
+    print(f"4. Custom GARCH(1,1) Jump Diffusion ({N:,} paths) [{DEVICE}]")
+    model = mbt.StochasticModel(
+        name="my_garch_jd",
+        drift="mu",
+        diffusion="sqrt(h)",
+        jump_intensity="lambda",
+        jump_size="normal(mu_j, sigma_j)",
+        state_vars={"h": 1e-4},
+        state_update={"h": "omega + alpha * (ret - mu) ** 2 + beta * h"},
+        params={
+            "mu": 0.08,
+            "omega": 1e-6,
+            "alpha": 0.10,
+            "beta": 0.85,
+            "lambda": 5.0,
+            "mu_j": -0.02,
+            "sigma_j": 0.04,
+        },
+    )
+    t0 = time.perf_counter()
+    result = mbt.run_stochastic(
+        model,
+        s0=100.0,
+        n_paths=N,
+        n_steps=252,
+        dt=1 / 252,
+        seed=42,
+        device=DEVICE,
+    )
+    elapsed = time.perf_counter() - t0
+    print(f"   Mean final price: {result['final_price']['mean']:.2f}")
+    print(f"   Max DD (P5):      {result['max_drawdown']['percentiles'][0][1]:.2%}")
+    print(f"   Elapsed: {elapsed:.3f}s\n")
+
+    # ── 5. Custom mean-reverting model (CPU — store_paths needs RAM) ────────
+    N_PLOT = 10_000
+    print(f"5. Custom mean-reverting model ({N_PLOT:,} paths) [cpu, store_paths]")
+    mean_rev = mbt.StochasticModel(
+        name="mean_reverting",
+        # Drift pulls price back toward 100
+        drift="kappa * (log(100.0) - log(S))",
+        diffusion="sigma",
+        params={"kappa": 2.0, "sigma": 0.25},
+    )
+    t0 = time.perf_counter()
+    result = mbt.run_stochastic(
+        mean_rev,
+        s0=80.0,   # start below mean
+        n_paths=N_PLOT,
+        n_steps=252,
+        dt=1 / 252,
+        seed=42,
+        store_paths=True,
+        device="cpu",
+    )
+    elapsed = time.perf_counter() - t0
+    print(f"   Mean final price: {result['final_price']['mean']:.2f} (target: 100)")
+    print(f"   Elapsed: {elapsed:.3f}s\n")
+
+    # ── 6. Fan chart visualization ──────────────────────────────────────────
+    print("6. Plotting fan chart...")
+    mbt.plot.stochastic_paths(
+        result,
+        title=f"Mean-reverting model (S0=80, target=100, {N_PLOT:,} paths)",
+        show=True,
+    )

examples/14_multi_timeframe.py

@@ -0,0 +1,84 @@
+"""Multi-Timeframe Strategy -- trend on 12h, entry on 1h.
+
+Demonstrates:
+  - bt.tf() for referencing higher-timeframe columns
+  - extra_timeframes config to inject resampled OHLCV
+  - Combining slow trend filter (12h EMA) with faster entry (1h RSI)
+
+Logic:
+  - 12h trend: EMA(20) > EMA(50) → bullish regime
+  - 1h entry:  RSI(14) < 35 during bullish regime → buy the dip
+  - Size: 50% of initial capital when conditions met, else flat
+
+Usage:
+    python examples/14_multi_timeframe.py
+"""
+import os
+import time
+import manifoldbt as mbt
+from manifoldbt.indicators import ema, rsi, close
+from manifoldbt.helpers import time_range, Slippage, Interval
+
+# -- Higher timeframe references ---------------------------------------------
+h12 = mbt.tf("12h")  # references columns like "12h.close"
+
+# -- Indicators ---------------------------------------------------------------
+# Trend filter on 12-hour bars (forward-filled onto 1h grid)
+trend_fast = ema(h12.close, 20)
+trend_slow = ema(h12.close, 50)
+bullish = trend_fast > trend_slow
+
+# Entry signal on 1-hour bars (native resolution)
+entry_rsi = rsi(close, 14)
+dip = entry_rsi < 35.0
+
+# -- Strategy -----------------------------------------------------------------
+strategy = (
+    mbt.Strategy.create("multi_tf_trend_dip")
+    .signal("bullish", bullish)
+    .signal("entry_rsi", entry_rsi)
+    .signal("dip", dip)
+    .size(mbt.when(mbt.col("bullish") & mbt.col("dip"), 0.5, 0.0))
+    .stop_loss(pct=3.0)
+    .describe("12h EMA trend + 1h RSI dip-buy, 3% stop-loss")
+)
+
+# -- Config -------------------------------------------------------------------
+start, end = time_range("2022-01-01", "2025-01-01")
+
+config = mbt.BacktestConfig(
+    universe=[1],
+    time_range_start=start,
+    time_range_end=end,
+    bar_interval=Interval.hours(1),
+    initial_capital=10_000,
+    execution=mbt.ExecutionConfig(
+        allow_short=False,
+        max_position_pct=0.5,
+        position_sizing_mode="FractionOfInitialCapital",
+    ),
+    fees=mbt.FeeConfig.binance_perps(),
+    slippage=Slippage.fixed_bps(2),
+    warmup_bars=50,
+    extra_timeframes={
+        "12h": Interval.hours(12),
+    },
+)
+
+# -- Run ----------------------------------------------------------------------
+if __name__ == "__main__":
+    root = os.path.join(os.path.dirname(__file__), "..")
+    store = mbt.DataStore(
+        data_root=os.path.abspath(os.path.join(root, "data")),
+        metadata_db=os.path.abspath(os.path.join(root, "metadata", "metadata.sqlite")),
+    )
+
+    t0 = time.perf_counter()
+    result = mbt.run(strategy, config, store)
+    elapsed = time.perf_counter() - t0
+    
+
+    print(result.summary())
+    print(f"\nElapsed: {elapsed:.3f}s")
+
+    mbt.plot.equity(result, show=True)

pyproject.toml

@@ -1,6 +1,6 @@
 [project]
 name = "manifoldbt"
-version = "0.2.0"
+version = "0.3.0"
 description = "Rust-powered backtesting engine for quantitative research"
 requires-python = ">=3.9"
 license = "MIT"

python/manifoldbt/__init__.py

@@ -24,6 +24,7 @@ from manifoldbt._native import (
     py_run_stability as _run_stability_native,
     py_replay as _replay_native,
     py_run_monte_carlo,
+    py_run_stochastic as _run_stochastic_native,
     run_portfolio as _run_portfolio_native,
     py_ingest as _ingest_native,
 )
@@ -42,7 +43,7 @@ from manifoldbt.exceptions import (
     LicenseError,
     StrategyError,
 )
-from manifoldbt.expr import AssetRef, Expr, asset, col, hold, lit, param, s, scan, symbol_ref, when
+from manifoldbt.expr import AssetRef, Expr, TimeframeRef, asset, col, hold, lit, param, s, scan, symbol_ref, tf, when
 from manifoldbt.helpers import (
     ExecutionPrice,
     FillModel,
@@ -75,9 +76,9 @@ def _print_banner():
         if tier == "Pro" and email:
             print(f"manifoldbt v{__version__} | \033[38;5;214mPro\033[0m | {email}")
         else:
-            print(f"manifoldbt v{__version__} | \033[36mCommunity\033[0m | upgrade: manifold-bt.com")
+            print(f"manifoldbt v{__version__} | \033[36mCommunity\033[0m | upgrade: www.manifoldbt.com")
     except Exception:
-        print(f"manifoldbt v{__version__} | \033[36mCommunity\033[0m | upgrade: manifold-bt.com")
+        print(f"manifoldbt v{__version__} | \033[36mCommunity\033[0m | upgrade: www.manifoldbt.com")
 
 _print_banner()
 del _print_banner
@@ -102,7 +103,7 @@ def _print_pro_summary() -> None:
         print()
         for w in _pro_warnings:
             print(f"\033[38;5;214m[!] {w} -- Pro feature\033[0m")
-        print("\033[38;5;214m  -> upgrade at manifold-bt.com\033[0m")
+        print("\033[38;5;214m  -> upgrade at www.manifoldbt.com\033[0m")
 
 
 import atexit
@@ -646,6 +647,105 @@ def replay(
     return Result(raw)
 
 
+# ---------------------------------------------------------------------------
+# Stochastic simulation API
+# ---------------------------------------------------------------------------
+
+from manifoldbt.stochastic import StochasticModel
+
+
+def run_stochastic(
+    model,
+    *,
+    s0: float = 100.0,
+    n_paths: int = 1000,
+    n_steps: int = 252,
+    dt: float = 1.0 / 252.0,
+    params: Optional[Dict[str, float]] = None,
+    seed: Optional[int] = None,
+    confidence_levels: Optional[List[float]] = None,
+    store_paths: bool = False,
+    device: str = "cpu",
+    precision: str = "f64",
+) -> Dict[str, Any]:
+    """Run a stochastic simulation via SDE expression DSL.
+
+    All expressions are compiled to native Rust and executed with Rayon
+    parallelism — no Python callback overhead.
+
+    Args:
+        model: Either a preset name (``"gbm"``, ``"heston"``, ``"merton"``,
+            ``"garch_jd"``) or a :class:`StochasticModel` instance.
+        s0: Initial price.
+        n_paths: Number of simulation paths.
+        n_steps: Number of time steps per path.
+        dt: Time step in years (``1/252`` = daily, ``1/252/390`` = minute).
+        params: Parameter overrides (merged with model defaults).
+        seed: RNG seed for reproducibility.
+        confidence_levels: Quantile levels for reporting.
+        store_paths: Whether to store full price paths.
+        device: ``"cpu"`` (default, Rayon parallel) or ``"cuda"``/``"gpu"``
+            (CUDA GPU, requires build with ``--features cuda``).
+        precision: ``"f64"`` (default, double) or ``"f32"`` (float, ~10-20x
+            faster on consumer GPUs, suitable for research/prototyping).
+
+    Returns:
+        Dict with ``final_price``, ``final_return``, ``max_drawdown``,
+        ``annualized_return``, ``annualized_vol`` (each with percentiles,
+        mean, std, min, max), and optionally ``paths`` (Arrow array) +
+        ``paths_n_steps``.
+
+    Example:
+        >>> result = mbt.run_stochastic("gbm", s0=100, n_paths=10000,
+        ...     n_steps=252, dt=1/252, params={"mu": 0.05, "sigma": 0.2})
+        >>> result["final_price"]["mean"]
+        105.12
+
+        >>> model = mbt.StochasticModel(
+        ...     drift="mu", diffusion="sqrt(h)",
+        ...     state_vars={"h": 1e-4},
+        ...     state_update={"h": "omega + alpha * (ret - mu)**2 + beta * h"},
+        ...     params={"mu": 0.08, "omega": 1e-6, "alpha": 0.1, "beta": 0.85},
+        ... )
+        >>> result = mbt.run_stochastic(model, s0=100, n_paths=5000)
+    """
+    config: Dict[str, Any] = {
+        "s0": s0,
+        "n_paths": n_paths,
+        "n_steps": n_steps,
+        "dt": dt,
+        "store_paths": store_paths,
+        "device": device,
+        "precision": precision,
+    }
+
+    if seed is not None:
+        config["rng_seed"] = seed
+
+    if confidence_levels is not None:
+        config["confidence_levels"] = confidence_levels
+
+    if isinstance(model, str):
+        # Preset name
+        config["preset"] = model
+        if params:
+            config["params"] = params
+    elif isinstance(model, StochasticModel):
+        model_dict = model.to_dict()
+        if params:
+            model_dict["params"].update(params)
+        config["model"] = model_dict
+    else:
+        raise TypeError(
+            f"model must be a preset name (str) or StochasticModel, got {type(model).__name__}"
+        )
+
+    try:
+        return _run_stochastic_native(json.dumps(config))
+    except (ValueError, RuntimeError) as exc:
+        raise _classify_error(exc) from exc
+
+
 # ---------------------------------------------------------------------------
 # Portfolio API
 # ---------------------------------------------------------------------------
@@ -748,6 +848,7 @@ __all__ = [
     # DSL
     "AssetRef",
     "Expr",
+    "TimeframeRef",
     "asset",
     "col",
     "lit",
@@ -755,6 +856,7 @@ __all__ = [
     "s",
     "scan",
     "symbol_ref",
+    "tf",
     "when",
     # Strategy & config
     "Strategy",
@@ -780,6 +882,9 @@ __all__ = [
     "run_stability",
     "replay",
     "py_run_monte_carlo",
+    # Stochastic simulation
+    "run_stochastic",
+    "StochasticModel",
     # Portfolio
     "Portfolio",
     "run_portfolio",

python/manifoldbt/_native.pyi

@@ -121,3 +121,6 @@ def py_run_monte_carlo(
     result: BacktestResult,
     mc_config_json: str,
 ) -> Dict[str, Any]: ...
+def py_run_stochastic(
+    sim_config_json: str,
+) -> Dict[str, Any]: ...

python/manifoldbt/config.py

@@ -191,6 +191,11 @@ class BacktestConfig:
     """When True, simulation runs on 1-minute bars regardless of bar_interval.
     Signals are still evaluated at bar_interval resolution (hybrid mode).
     Use for precise SL/TP fills and intraday drawdown tracking. Slower."""
+    extra_timeframes: Dict[str, Any] = field(default_factory=dict)
+    """Additional timeframes for multi-timeframe strategies.
+    Maps labels to Interval dicts. The engine resamples native bars
+    and injects prefixed columns (e.g. "1h.close", "4h.high").
+    Example: ``{"1h": Interval.hours(1), "4h": Interval.hours(4)}``"""
 
     def to_json_dict(self) -> dict:
         d: dict = {
@@ -219,6 +224,8 @@ class BacktestConfig:
             d["symbol_names"] = self.symbol_names
         if self.warmup_bars > 0:
             d["warmup_bars"] = self.warmup_bars
+        if self.extra_timeframes:
+            d["extra_timeframes"] = self.extra_timeframes
         return d
 
     def to_json(self) -> str:

python/manifoldbt/expr.py

@@ -560,6 +560,62 @@ def asset(symbol: str) -> AssetRef:
     return AssetRef(symbol)
 
 
+class TimeframeRef:
+    """Reference columns from a higher timeframe.
+
+    The columns are forward-filled: a completed 1h bar's value becomes
+    available at the start of the *next* 1h bar and persists until that
+    bar completes.  This avoids lookahead bias.
+
+    Requires ``extra_timeframes`` in ``BacktestConfig``.
+    """
+
+    __slots__ = ("_tf",)
+
+    def __init__(self, tf: str) -> None:
+        self._tf = tf
+
+    @property
+    def open(self) -> Expr:
+        return col(f"{self._tf}.open")
+
+    @property
+    def high(self) -> Expr:
+        return col(f"{self._tf}.high")
+
+    @property
+    def low(self) -> Expr:
+        return col(f"{self._tf}.low")
+
+    @property
+    def close(self) -> Expr:
+        return col(f"{self._tf}.close")
+
+    @property
+    def volume(self) -> Expr:
+        return col(f"{self._tf}.volume")
+
+    def col(self, name: str) -> Expr:
+        """Reference any column from this timeframe."""
+        return col(f"{self._tf}.{name}")
+
+    def __repr__(self) -> str:
+        return f"TimeframeRef({self._tf!r})"
+
+
+def tf(timeframe: str) -> TimeframeRef:
+    """Reference a higher timeframe for multi-TF strategies.
+
+    Usage::
+
+        h1 = bt.tf("1h")
+        trend = ema(h1.close, 20) > ema(h1.close, 50)
+
+    Requires ``extra_timeframes={"1h": Interval.hours(1)}`` in config.
+    """
+    return TimeframeRef(timeframe)
+
+
 # ---------------------------------------------------------------------------
 # Scan (stateful fold) support
 # ---------------------------------------------------------------------------

python/manifoldbt/plot/__init__.py

@@ -43,6 +43,7 @@ from manifoldbt.plot.research import (
     heatmap_2d,
     monte_carlo,
     stability,
+    stochastic_paths,
     surface_3d,
     walk_forward,
 )
@@ -73,6 +74,7 @@ __all__ = [
     "stability",
     "correlation_matrix",
     "monte_carlo",
+    "stochastic_paths",
     # Composites
     "tearsheet",
     "research_report",

python/manifoldbt/plot/research.py

@@ -717,3 +717,122 @@ def monte_carlo(
         ax_.set_ylabel("Equity")
         ax_.legend(loc="upper left", fontsize=7, framealpha=0.3)
         return finalize(fig, show=show, save=save)
+
+
+# ── Stochastic Simulation Paths ───────────────────────────────────────────
+
+
+def stochastic_paths(
+    result: Dict[str, Any],
+    *,
+    percentiles: Optional[List[int]] = None,
+    n_sample_paths: int = 50,
+    ax: Optional[Axes] = None,
+    median_color: str = ACCENT,
+    band_color: str = ACCENT,
+    title: Optional[str] = None,
+    figsize: Tuple[float, float] = (12, 5),
+    show: bool = False,
+    save: Optional[Union[str, Path]] = None,
+) -> Figure:
+    """Fan chart for stochastic simulation paths with percentile bands.
+
+    Args:
+        result: Dict returned by ``mbt.run_stochastic(..., store_paths=True)``.
+            Must contain ``paths`` (flat Arrow array) and ``paths_n_steps``.
+        percentiles: Percentile levels for bands. Default ``[5, 25, 50, 75, 95]``.
+        n_sample_paths: Number of individual paths to draw (faded). 0 to disable.
+    """
+    if percentiles is None:
+        percentiles = [5, 25, 50, 75, 95]
+
+    paths_raw = result.get("paths")
+    n_steps = result.get("paths_n_steps")
+    n_paths = result.get("n_paths", 0)
+    model_name = result.get("model_name", "stochastic")
+
+    if paths_raw is None or n_steps is None:
+        raise ValueError(
+            "result has no paths data. Run with store_paths=True."
+        )
+
+    # Reshape flat Arrow/numpy array → (n_paths, n_steps+1)
+    flat = np.asarray(paths_raw, dtype=np.float64)
+    paths = flat.reshape((n_paths, n_steps))
+
+    if title is None:
+        title = f"Stochastic simulation - {model_name} ({n_paths:,} paths)"
+
+    with theme_context():
+        fig, ax_ = get_or_create_ax(ax, figsize)
+
+        x = np.arange(paths.shape[1])
+
+        # Draw sample paths (faded)
+        if n_sample_paths > 0:
+            for i in range(min(n_sample_paths, n_paths)):
+                ax_.plot(x, paths[i], color=band_color, linewidth=0.3, alpha=0.06)
+
+        # Compute percentile bands
+        pct_lines = {pct: np.percentile(paths, pct, axis=0) for pct in percentiles}
+
+        # Fill between symmetric bands
+        for lo, hi in [(0, -1), (1, -2)]:
+            if lo < len(percentiles) and hi < 0 and abs(hi) <= len(percentiles):
+                ax_.fill_between(
+                    x,
+                    pct_lines[percentiles[lo]],
+                    pct_lines[percentiles[hi]],
+                    color=band_color,
+                    alpha=0.08,
+                )
+
+        # Percentile lines
+        s0 = paths[0, 0] if paths.shape[1] > 0 else 100.0
+        for pct in percentiles:
+            final = pct_lines[pct][-1]
+            ret_pct = (final / s0 - 1) * 100
+            if pct == 50:
+                ax_.plot(
+                    x, pct_lines[pct], color=median_color, linewidth=2,
+                    label=f"P{pct} (median): {ret_pct:+.1f}%", zorder=3,
+                )
+            else:
+                ax_.plot(
+                    x, pct_lines[pct], color=band_color, linewidth=0.5,
+                    alpha=0.4, label=f"P{pct}: {ret_pct:+.1f}%",
+                )
+
+        # Stats box
+        final_prices = paths[:, -1]
+        running_peak = np.maximum.accumulate(paths, axis=1)
+        drawdowns = (paths - running_peak) / running_peak
+        max_dd_per_path = drawdowns.min(axis=1) * 100
+
+        dd_p5 = np.percentile(max_dd_per_path, 5)
+        dd_p50 = np.percentile(max_dd_per_path, 50)
+        mean_ret = (np.mean(final_prices) / s0 - 1) * 100
+
+        stats_text = (
+            f"Mean return: {mean_ret:+.1f}%\n"
+            f"Max DD (P5): {dd_p5:.1f}%\n"
+            f"Max DD (P50): {dd_p50:.1f}%"
+        )
+        ax_.text(
+            0.98, 0.95, stats_text,
+            transform=ax_.transAxes, ha="right", va="top",
+            color="#8a8a8a", fontsize=8, fontfamily="monospace",
+            bbox={
+                "boxstyle": "round,pad=0.4",
+                "facecolor": "#111116",
+                "edgecolor": "#1e1e24",
+                "alpha": 0.9,
+            },
+        )
+
+        ax_.margins(x=0.02)
+        ax_.set_title(title)
+        ax_.set_xlabel("Time steps")
+        ax_.set_ylabel("Price")
+        ax_.legend(loc="upper left", fontsize=7, framealpha=0.3)
+        return finalize(fig, show=show, save=save)

python/manifoldbt/stochastic.py

@@ -0,0 +1,177 @@
+"""Stochastic simulation via SDE expression DSL.
+
+Define stochastic differential equations as string expressions and simulate
+price paths at native Rust speed with Rayon parallelism.
+
+Example
+-------
+>>> import manifoldbt as mbt
+>>> model = mbt.StochasticModel(
+...     drift="mu",
+...     diffusion="sqrt(h)",
+...     jump_intensity="lambda",
+...     jump_size="normal(mu_j, sigma_j)",
+...     state_vars={"h": 1e-4},
+...     state_update={"h": "omega + alpha * (ret - mu) ** 2 + beta * h"},
+...     params={"mu": 0.08, "omega": 1e-6, "alpha": 0.1, "beta": 0.85,
+...             "lambda": 5.0, "mu_j": -0.02, "sigma_j": 0.04},
+... )
+>>> result = mbt.run_stochastic(model, s0=100, n_paths=10000, n_steps=252, dt=1/252)
+
+Presets
+-------
+>>> result = mbt.run_stochastic("gbm", s0=100, n_paths=10000, n_steps=252, dt=1/252,
+...                              params={"mu": 0.05, "sigma": 0.2})
+"""
+
+from typing import Any, Dict, List, Optional
+
+
+class StochasticModel:
+    """Define a stochastic differential equation model via string expressions.
+
+    The SDE has the form::
+
+        dS = drift(S,t,state) * S * dt
+           + diffusion(S,t,state) * S * dW
+           + jump_size * dN(jump_intensity)
+
+    where state variables (like GARCH variance ``h``) are updated after each
+    price step according to ``state_update`` expressions.
+
+    Parameters
+    ----------
+    drift : str
+        Drift expression (μ). E.g. ``"mu"`` or ``"mu - 0.5 * h"``.
+    diffusion : str
+        Diffusion expression (σ). E.g. ``"sigma"`` or ``"sqrt(h)"``.
+    jump_intensity : str, optional
+        Jump intensity expression (λ). E.g. ``"lambda"``.
+    jump_size : str, optional
+        Jump size expression. E.g. ``"normal(mu_j, sigma_j)"``.
+    state_vars : dict, optional
+        State variable names and initial values. ``S`` is implicit (index 0).
+    state_update : dict, optional
+        Update expressions for state variables.
+        E.g. ``{"h": "omega + alpha * (ret - mu) ** 2 + beta * h"}``.
+    params : dict, optional
+        Model parameter names and values.
+    name : str
+        Model name (default ``"custom"``).
+
+    Available identifiers in expressions
+    -------------------------------------
+    - Any key in ``params`` → parameter value
+    - Any key in ``state_vars`` → current state value
+    - ``S`` → current price
+    - ``ret`` → log-return from previous step
+    - ``dt`` → time step size
+    - ``t`` → current simulation time
+    - ``step`` → current step index
+
+    Available functions
+    -------------------
+    ``sqrt``, ``abs``, ``log``, ``exp``, ``floor``, ``max``, ``min``, ``pow``,
+    ``normal(mu, sigma)``, ``uniform(lo, hi)``, ``randn()``,
+    ``if(cond, then, else)``
+
+    Operators: ``+``, ``-``, ``*``, ``/``, ``**``, ``>``, ``<``, ``>=``,
+    ``<=``, ``==``, ``&&``, ``||``, ``!``, ternary ``cond ? a : b``
+    """
+
+    def __init__(
+        self,
+        *,
+        drift: str,
+        diffusion: str,
+        jump_intensity: Optional[str] = None,
+        jump_size: Optional[str] = None,
+        state_vars: Optional[Dict[str, float]] = None,
+        state_update: Optional[Dict[str, str]] = None,
+        params: Optional[Dict[str, float]] = None,
+        name: str = "custom",
+    ):
+        self.name = name
+        self.drift = drift
+        self.diffusion = diffusion
+        self.jump_intensity = jump_intensity
+        self.jump_size = jump_size
+        self.state_vars = state_vars or {}
+        self.state_update = state_update or {}
+        self.params = params or {}
+
+    def to_dict(self) -> Dict[str, Any]:
+        """Serialise to a dict suitable for JSON encoding."""
+        d: Dict[str, Any] = {
+            "name": self.name,
+            "drift": self.drift,
+            "diffusion": self.diffusion,
+            "params": dict(self.params),
+        }
+        if self.jump_intensity is not None:
+            d["jump_intensity"] = self.jump_intensity
+        if self.jump_size is not None:
+            d["jump_size"] = self.jump_size
+        if self.state_vars:
+            d["state_vars"] = dict(self.state_vars)
+        if self.state_update:
+            d["state_update"] = dict(self.state_update)
+        return d
+
+    def __repr__(self) -> str:
+        parts = [f"StochasticModel(name={self.name!r}"]
+        parts.append(f"drift={self.drift!r}")
+        parts.append(f"diffusion={self.diffusion!r}")
+        if self.jump_intensity:
+            parts.append(f"jump_intensity={self.jump_intensity!r}")
+        if self.state_vars:
+            parts.append(f"state_vars={self.state_vars!r}")
+        return ", ".join(parts) + ")"
+
+
+# Preset shortcuts
+GBM = StochasticModel(
+    name="gbm",
+    drift="mu",
+    diffusion="sigma",
+    params={"mu": 0.05, "sigma": 0.2},
+)
+
+HESTON = StochasticModel(
+    name="heston",
+    drift="mu",
+    diffusion="sqrt(v)",
+    state_vars={"v": 0.04},
+    state_update={
+        "v": "max(v + kappa * (theta - v) * dt + xi * sqrt(v) * sqrt(dt) * randn(), 0.0)"
+    },
+    params={"mu": 0.05, "kappa": 2.0, "theta": 0.04, "xi": 0.3},
+)
+
+MERTON = StochasticModel(
+    name="merton",
+    drift="mu",
+    diffusion="sigma",
+    jump_intensity="lambda",
+    jump_size="normal(mu_j, sigma_j)",
+    params={"mu": 0.05, "sigma": 0.2, "lambda": 1.0, "mu_j": -0.05, "sigma_j": 0.08},
+)
+
+GARCH_JD = StochasticModel(
+    name="garch_jd",
+    drift="mu",
+    diffusion="sqrt(h)",
+    jump_intensity="lambda",
+    jump_size="normal(mu_j, sigma_j)",
+    state_vars={"h": 1e-4},
+    state_update={"h": "omega + alpha * (ret - mu) ** 2 + beta * h"},
+    params={
+        "mu": 0.08,
+        "omega": 1e-6,
+        "alpha": 0.1,
+        "beta": 0.85,
+        "lambda": 5.0,
+        "mu_j": -0.02,
+        "sigma_j": 0.04,
+    },
+)