Simulation & Optimization

Monte Carlo portfolio simulation, Bayesian parameter optimization, and stress testing. Rust-native PRNG and GP surrogate.

Pro Feature. Requires a Pro or Ultra subscription. Get started at api.mathematicalcompany.com

Monte Carlo Simulation

Horizon includes a Monte Carlo simulation engine for stress-testing prediction market portfolios. The entire simulation runs in Rust with a custom PRNG (xoshiro256++), Box-Muller normal generation, and Cholesky decomposition for correlated outcomes. No external dependencies.

Monte Carlo simulation answers the question: “Given my current positions and probability estimates, what is the distribution of possible portfolio outcomes?” This helps you understand tail risk (VaR, CVaR), win probability, and the impact of correlation between positions.

Overview

Fast Simulation

50,000 scenarios in milliseconds. All math runs in Rust with zero Python overhead.

Correlated Outcomes

Model dependencies between markets via correlation matrices. Uses Cholesky decomposition.

Full Risk Stats

VaR (95/99), CVaR, win probability, max loss/gain, percentile distribution, and every scenario PnL.

Deterministic Seeds

Reproducible results with seed parameter. Same seed = identical output every time.

---

Quick Start

import horizon as hz

# Define your positions
positions = [
    hz.SimPosition("election-winner", "yes", 100.0, 0.50, 0.60),
    hz.SimPosition("fed-rate-cut", "no", 50.0, 0.40, 0.30),
]

# Run simulation
result = hz.simulate(positions=positions, scenarios=50000, seed=42)

print(f"Mean PnL:    ${result.mean_pnl:.2f}")
print(f"VaR 95:      ${result.var_95:.2f}")
print(f"CVaR 95:     ${result.cvar_95:.2f}")
print(f"Win Prob:    {result.win_probability:.1%}")
print(f"Max Loss:    ${result.max_loss:.2f}")
print(f"Max Gain:    ${result.max_gain:.2f}")

---

SimPosition

Represents a position for simulation input.

pos = hz.SimPosition(
    market_id="election-winner",  # Market identifier
    side="yes",                    # "yes" or "no" (case insensitive)
    size=100.0,                    # Position size in contracts
    entry_price=0.50,              # Average entry price
    current_price=0.60,            # Current market probability estimate
)

Parameter	Type	Description
market_id	str	Market identifier
side	str	"yes" or "no" (case insensitive)
size	float	Position size in contracts
entry_price	float	Average entry price
current_price	float	Current market probability / price

The current_price is used as the true probability for the Bernoulli draw:

• Yes position: wins size * (1 - entry_price) if outcome = 1, loses size * entry_price if outcome = 0
• No position: wins size * entry_price if outcome = 0, loses size * (1 - entry_price) if outcome = 1

SimulationResult

The result object contains full risk statistics.

result = hz.monte_carlo(positions, scenarios=10000, seed=42)

Field	Type	Description
mean_pnl	float	Average PnL across all scenarios
median_pnl	float	Median PnL
std_dev	float	Standard deviation of PnL
var_95	float	Value at Risk (5th percentile)
var_99	float	Value at Risk (1st percentile)
cvar_95	float	Conditional VaR (expected shortfall at 5%)
max_loss	float	Worst-case scenario PnL
max_gain	float	Best-case scenario PnL
win_probability	float	Fraction of scenarios with positive PnL
scenario_pnl	list[float]	All scenario PnLs (sorted ascending)
percentiles	list[tuple[float, float]]	Percentile distribution (percentile, value)

---

Core Functions

hz.monte_carlo

Run the Monte Carlo simulation directly (Rust function).

result = hz.monte_carlo(
    positions=positions,         # list[SimPosition]
    scenarios=10000,             # Number of simulation runs
    correlation_matrix=None,     # Optional NxN correlation matrix
    seed=42,                     # Random seed (None = random)
)

hz.simulate

Python wrapper with ergonomic features.

result = hz.simulate(
    engine=None,                 # Engine (auto-extract positions)
    positions=None,              # Explicit positions (overrides engine)
    scenarios=10000,             # Number of scenarios
    correlations=None,           # dict[tuple[str,str], float] correlation pairs
    prices=None,                 # dict[str, float] override current prices
    seed=None,                   # Random seed
)

The wrapper provides two key conveniences:

Auto-extraction from Engine: If engine is passed, positions are extracted from engine.positions() and current prices from engine.all_feed_snapshots().

Correlation dict: Instead of building a full NxN matrix, pass a dict of pairwise correlations:

result = hz.simulate(
    positions=positions,
    correlations={
        ("election-winner", "fed-rate-cut"): 0.3,
        ("election-winner", "recession"): -0.5,
    },
)

The wrapper builds the full correlation matrix automatically, including handling reverse key order (both ("a", "b") and ("b", "a") work).

---

Examples

Portfolio Risk Assessment

import horizon as hz

positions = [
    hz.SimPosition("trump-win", "yes", 200.0, 0.45, 0.55),
    hz.SimPosition("fed-cut-march", "yes", 100.0, 0.60, 0.65),
    hz.SimPosition("btc-100k", "no", 150.0, 0.30, 0.25),
]

result = hz.simulate(
    positions=positions,
    correlations={
        ("trump-win", "fed-cut-march"): 0.2,
        ("trump-win", "btc-100k"): -0.1,
    },
    scenarios=50000,
    seed=42,
)

print(f"Expected PnL:  ${result.mean_pnl:+.2f}")
print(f"Std Dev:        ${result.std_dev:.2f}")
print(f"VaR 95:         ${result.var_95:+.2f}")
print(f"CVaR 95:        ${result.cvar_95:+.2f}")
print(f"Win Rate:       {result.win_probability:.1%}")
print(f"Best Case:      ${result.max_gain:+.2f}")
print(f"Worst Case:     ${result.max_loss:+.2f}")

Simulate from Live Engine

engine = hz.Engine(risk_config=hz.RiskConfig(max_position_per_market=200))

# ... trade for a while, build up positions ...

# Stress test current portfolio
result = hz.simulate(engine=engine, scenarios=50000)
print(f"Portfolio VaR 95: ${result.var_95:.2f}")

Comparing Correlated vs Uncorrelated Risk

import horizon as hz

positions = [
    hz.SimPosition("mkt1", "yes", 100.0, 0.50, 0.60),
    hz.SimPosition("mkt2", "yes", 100.0, 0.50, 0.60),
]

# Uncorrelated
r_uncorr = hz.simulate(positions=positions, scenarios=50000, seed=42)

# Highly correlated
r_corr = hz.simulate(
    positions=positions,
    correlations={("mkt1", "mkt2"): 0.9},
    scenarios=50000,
    seed=42,
)

print(f"Uncorrelated StdDev: ${r_uncorr.std_dev:.2f}")
print(f"Correlated StdDev:   ${r_corr.std_dev:.2f}")
# Correlated positions have higher portfolio variance

Seed Determinism

# Same seed = identical results
r1 = hz.monte_carlo(positions, 10000, None, 42)
r2 = hz.monte_carlo(positions, 10000, None, 42)
assert r1.mean_pnl == r2.mean_pnl  # Exact match

# Different seed = different (but statistically similar) results
r3 = hz.monte_carlo(positions, 10000, None, 99)
# r3.mean_pnl will be close to r1.mean_pnl but not identical

---

Mathematical Background

Binary Outcome PnL

Each position is a binary contract. If the event resolves YES (outcome = 1):

• YES position PnL = size * (1 - entry_price)
• NO position PnL = size * (-(1 - entry_price))

If the event resolves NO (outcome = 0):

• YES position PnL = size * (-entry_price)
• NO position PnL = size * entry_price

The current_price is the probability of outcome = 1 in each simulation draw.

Correlated Draws via Cholesky

To model correlated binary outcomes:

1. Draw N standard normal variables (Box-Muller transform)
2. Apply Cholesky decomposition of the correlation matrix to induce correlations
3. Convert correlated normals to uniform via the normal CDF
4. Threshold each uniform against the position’s current_price to produce a Bernoulli outcome

This produces correlated binary outcomes that respect the specified correlation structure while maintaining the correct marginal probabilities.

VaR and CVaR

• VaR 95: The 5th percentile of the PnL distribution. “With 95% confidence, losses will not exceed this amount.”
• VaR 99: The 1st percentile. More conservative.
• CVaR 95 (Expected Shortfall): The average PnL of the worst 5% of scenarios. Always ≤ VaR 95.

Bayesian Optimization

Find optimal strategy parameters using a Gaussian Process surrogate with Expected Improvement acquisition.

from horizon import bayesian_optimize

# Define a pipeline factory: takes params, returns pipeline list
def make_pipeline(params):
    return [
        hz.market_maker(feed_name="book", gamma=params["gamma"], size=params["size"]),
    ]

result = bayesian_optimize(
    data=historical_data,          # same format as hz.backtest()
    pipeline_factory=make_pipeline,
    search_space={
        "gamma": (0.1, 1.0),
        "size": (1.0, 20.0),
    },
    objective="sharpe_ratio",      # metric to maximize from BacktestResult
    n_trials=50,
    n_initial=10,
    seed=42,
)

print(result.best_params)   # {"gamma": 0.5, "size": 8.2}
print(result.best_score)    # 1.85 (best Sharpe found)
print(result.n_trials)      # 50
print(result.all_trials)    # List of {"params": ..., "score": ...}
print(result.convergence)   # Best score at each trial index

Zero external dependencies. Uses RBF kernel GP with Cholesky decomposition.

Monte Carlo simulation assumes current_price represents the true probability. If your price estimates are poor, the simulation results will be misleading. Consider running simulations across a range of probability assumptions to understand sensitivity.