var-mc — Monte Carlo Value at Risk in Rust

A fast, correct, reproducible Monte Carlo Value at Risk (VaR) engine written in Rust with optional Python bindings. It simulates correlated multi-asset portfolio returns and computes VaR, Conditional VaR (Expected Shortfall), and an analytic risk-contribution breakdown (component & marginal VaR) — a few times faster than well-optimized NumPy, with true multi-core scaling and deterministic seeding.

Design goals: get the risk numbers right, make them reproducible, and make them fast — in that order. The library is validated against closed-form Gaussian VaR and against a NumPy reference implementation, and the tricky parts (covariance scaling, thread-independent seeding) are pinned by regression tests.

---

Table of contents

• What it computes
• Quick start (Python)
• Quick start (Rust)
• Methodology
• Correctness & testing
• Performance
• Feature matrix
• Building from source
• Project layout
• References
• License

---

What it computes

Given a portfolio value, asset weights, annualized mean returns, and an annualized covariance matrix, the engine simulates the portfolio's return distribution over a chosen horizon and reports:

Metric	Meaning
VaR at 95% / 99% (or any level)	the loss that is exceeded only 1 − α of the time
CVaR / Expected Shortfall	the average loss conditional on being in the tail beyond VaR
Parametric VaR	closed-form delta-normal VaR, reported alongside the MC estimate as a sanity cross-check
Component VaR	each asset's additive contribution to portfolio VaR (sums to the parametric VaR)
Marginal VaR	sensitivity of portfolio VaR to each asset's weight, ∂VaR/∂wᵢ

Two innovation distributions are supported — Gaussian and standardized Student-t (heavier tails) — plus optional antithetic variates for variance reduction.

---

Quick start (Python)

pip install maturin numpy
maturin develop --release      # builds the Rust extension into your virtualenv

import numpy as np
import var_mc

# Annualized inputs for a 3-asset portfolio
weights = np.array([0.5, 0.3, 0.2])
returns = np.array([0.08, 0.12, 0.05])          # expected annual returns
cov     = np.array([[0.0400, 0.0180, 0.0030],   # annual covariance
                    [0.0180, 0.0900, 0.0090],
                    [0.0030, 0.0090, 0.0225]])

# 1) One-shot functional API
res = var_mc.calculate_var(
    portfolio_value=10_000_000,
    weights=weights, returns=returns, covariance=cov,
    time_horizon=10,                 # trading days
    confidence_levels=[0.95, 0.99],
    n_simulations=1_000_000,
)
print(f"10-day 95% VaR:  ${res['var_95']:,.0f}")
print(f"10-day 99% CVaR: ${res['cvar_99']:,.0f}")
# res also has: parametric_var_95/99, execution_time, memory_usage (estimated MB)

# 2) Object API (standard 95%/99%)
calc = var_mc.VaRCalculator()
r = calc.calculate(portfolio_value=10_000_000, weights=weights,
                   returns_mean=returns, returns_cov=cov,
                   time_horizon=10, n_simulations=1_000_000)
print(r)   # VaRResult(var_95=..., var_99=..., cvar_95=..., cvar_99=...)

# 3) Reusable engine: heavy tails + variance reduction + risk decomposition
engine = var_mc.MonteCarloEngine(
    n_simulations=1_000_000, time_horizon=10,
    distribution="t", degrees_of_freedom=5,   # Student-t innovations
    antithetic=True,
)
tail = engine.calculate_var(10_000_000, weights, returns, cov,
                            confidence_levels=[0.95, 0.99])
contrib = engine.risk_contributions(10_000_000, weights, returns, cov,
                                    confidence_level=0.95)
print("Component VaR:", contrib["component_var"])   # sums to parametric VaR
print("Marginal VaR:", contrib["marginal_var"])

calculate_var returns a dict keyed by integer confidence percent (var_95, cvar_95, parametric_var_95, ...) plus execution_time (seconds) and memory_usage (estimated simulation-buffer size in MB).

---

Quick start (Rust)

# Cargo.toml
[dependencies]
var-mc = { git = "https://github.com/yourorg/var-monte-carlo" }
ndarray = "0.15"

use ndarray::array;
use var_mc::{run_var, Distribution, Portfolio, SimulationConfig};

let weights = array![0.5, 0.3, 0.2];
let returns = array![0.08, 0.12, 0.05];
let cov = array![
    [0.0400, 0.0180, 0.0030],
    [0.0180, 0.0900, 0.0090],
    [0.0030, 0.0090, 0.0225],
];

let portfolio = Portfolio::new(weights, returns, cov)?;   // validates dims, weights, PD
let config = SimulationConfig {
    n_simulations: 1_000_000,
    time_horizon: 10,
    seed: Some(42),
    parallel: true,
    distribution: Distribution::Normal,
    antithetic: false,
};
let report = run_var(config, &portfolio, 10_000_000.0, &[0.95, 0.99])?;
println!("95% VaR = {:.0}", report.var[0]);
println!("99% CVaR = {:.0}", report.cvar[1]);
# Ok::<(), var_mc::MonteCarloError>(())

The pure-Rust core has no Python dependency and can be used as a normal crate.

---

Methodology

For each of N simulations the engine draws a standard-normal vector z, correlates it with the Cholesky factor L of the daily covariance, scales it to the horizon, and adds the drift:

Σ_daily = Σ_annual / 252            μ_daily = μ_annual / 252
L        = chol(Σ_daily)            (lower-triangular, L·Lᵀ = Σ_daily)

rᵢ = μ_daily·T  +  (L·z)·√T·s       s = 1 (Normal)  or  √((df−2)/W) (Student-t, W~χ²_df)

Portfolio returns are r · w, losses are −V · (r·w), and VaR/CVaR are read off the sorted loss distribution. Component and marginal VaR are computed analytically (delta-normal), which is exact and — by construction — the component contributions sum to the parametric portfolio VaR.

Why "daily"? Annualized covariance must be scaled to the simulation step before the Cholesky, not after. Getting this wrong inflates VaR by ≈ √252 ≈ 16×. That exact mistake is guarded by a regression test (see below).

---

Correctness & testing

Numerical correctness is the whole point of a risk library, so it is tested first:

• MC vs. closed form — Gaussian Monte Carlo VaR matches the analytic V·(z_α·σ_T − μ_T) within Monte Carlo error at both 95% and 99%.
• Covariance scaling — the empirical volatility of simulated returns equals the horizon volatility (this is the √252 units-bug guard).
• Thread-independent reproducibility — the same seed produces identical results whether run on 1 core or 32; parallel and sequential paths are asserted equal.
• Heavy tails — Student-t VaR at 99% strictly exceeds Gaussian VaR.
• Risk decomposition — component VaR sums to the parametric VaR; marginal × weight equals component, element-wise.
• Variance reduction — with antithetic variates and an even sample count, the +z/−z shocks cancel exactly, so each asset's sample mean equals its drift.
• Input validation — non-conforming weights, dimension mismatches, and non-positive-definite covariances are rejected with typed errors.

cargo test                                   # pure-Rust core, no Python needed
cargo clippy --all-targets -- -D warnings    # lint clean

Cross-checking against NumPy (a properly Cholesky-based reference) agrees to within < 1% — the residual is ordinary Monte Carlo sampling error, not bias.

---

Performance

Measured on a 16-core machine, 1,000,000 simulations, 10-day horizon, release build (benchmark_test.py). The honest baseline is NumPy that uses the same Cholesky-based multi-asset model (numpy_optimized), which is already BLAS-vectorized:

Assets	NumPy (Cholesky)	var-mc (Rust)	Speedup	VaR agreement
10	0.151 s	0.031 s	~4.9×	< 0.3%
50	0.643 s	0.148 s	~4.3×	< 0.1%
100	1.269 s	0.247 s	~5.1×	< 0.2%
500	6.556 s	2.712 s	~2.4×	< 0.9%

Across the whole grid Rust averaged ~9 million scenarios/second. The advantage is largest for the realistic 100K–10M-simulation workloads; for tiny runs (≤ 10K sims) the fixed setup cost (validation, Cholesky, thread pool) isn't amortized and Rust can be on par with or slower than NumPy — Monte Carlo isn't the tool for 10,000 draws.

A note on honesty: you will sometimes see "50–100× faster than Python" claims for this kind of workload. Those hold only against pure-Python loops. Against vectorized, BLAS-backed NumPy doing the identical computation, a realistic figure is a few times faster — which is what you see above. The advantage grows for small-to-mid portfolios and RNG-bound workloads, and narrows for very large matrices where NumPy's OpenBLAS gemm dominates. The real wins beyond raw speed are multi-core scaling without the GIL, deterministic reproducibility, heavy-tailed simulation, and an embeddable zero-dependency Rust core.

The numpy_basic row in the benchmark models the whole portfolio as a single aggregated normal (O(sims), not O(sims × assets)) and is not an equivalent computation — it is included for context, not as a speedup baseline.

Reproduce it yourself:

pip install numpy pandas matplotlib   # numba optional
maturin develop --release
python benchmark_test.py

---

Feature matrix

Feature	Status
Correlated multi-asset Monte Carlo (Cholesky)	✅
VaR + CVaR / Expected Shortfall, arbitrary confidence levels	✅
Component & marginal VaR (delta-normal)	✅
Gaussian and Student-t innovations	✅
Antithetic variates (variance reduction)	✅
Multi-core parallelism (rayon), GIL released in Python	✅
Deterministic, thread-count-independent seeding	✅
Parametric VaR cross-check	✅
Quasi-Monte Carlo / Sobol sequences	🔜
Explicit SIMD kernels & memory pooling	🔜
Historical simulation / bootstrapping	🔜
Zero-copy NumPy output (results are currently copied)	🔜

---

Building from source

Requirements: a stable Rust toolchain, Python ≥ 3.9, and maturin.

git clone https://github.com/yourorg/var-monte-carlo
cd var-monte-carlo

# Rust core
cargo test
cargo build --release

# Python extension (into the active virtualenv)
python -m venv .venv && source .venv/bin/activate
pip install maturin numpy
maturin develop --release
python -c "import var_mc; print('ok')"

CI (GitHub Actions, see .github/workflows/ci.yml) runs cargo fmt, clippy, and cargo test, then builds and smoke-tests the wheel on Linux and macOS across Python 3.11–3.13.

---

Project layout

src/
├── lib.rs          crate root; re-exports; `python` bindings behind a cargo feature
├── portfolio.rs    Portfolio type + validation (weights, dims, positive-definiteness)
├── monte_carlo.rs  simulation engine: daily-cov Cholesky, chunked parallel RNG,
│                     Normal & Student-t, antithetic variates
├── var.rs          VaR / CVaR / component / marginal VaR + run_var
├── stats.rs        inverse normal CDF (Acklam) for the delta-normal metrics
└── python.rs       PyO3 + rust-numpy bindings (compiled with --features python)
tests/integration.rs   numerics & reproducibility regression tests
benchmark_test.py      NumPy / Numba / multiprocessing / Rust comparison suite

The PyO3/NumPy dependencies are optional (behind the python feature) so cargo test exercises the pure-Rust core without linking libpython.

---

References

1. Jorion, P. (2006). Value at Risk: The New Benchmark for Managing Financial Risk.
2. Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering.
3. McNeil, A. J., Frey, R., & Embrechts, P. (2015). Quantitative Risk Management.
4. Acklam, P. J. An algorithm for computing the inverse normal cumulative distribution function.

---

License

Dual-licensed under either of MIT or Apache-2.0, at your option.