@veenie/risk

A dependency-free Monte Carlo & Latin Hypercube risk simulation engine in TypeScript. 17 probability distributions, Iman–Conover rank correlation (the @RISK method), machine-precision quantile functions, and finance statistics (VaR, CVaR, tornado sensitivity) — in one small package with zero runtime dependencies.

Built to power MonteSheet, the Monte Carlo simulator for Google Sheets. Play with every distribution live at montesheet.com/distributions — interactive sandboxes running 10,000 simulations of this exact package in your browser.

Sibling package: @veenie/aero (lifting-line wing aerodynamics, powers AeroSheet). Both from veenie.space.

Quick start

npm install @veenie/risk

import { generateCorrelatedSamples, computeStats } from '@veenie/risk';

const variables = [
  { distId: 'normal', base: 100, params: { spread: 0.1 } },          // revenue
  { distId: 'pert',   base: 50,  params: { min: 30, mode: 50, max: 90 } } // cost
];

// Optional: how they move together
const correlation = [
  [1.0, 0.8],
  [0.8, 1.0]
];

// 100k iterations, reproducible, in milliseconds
const samples = generateCorrelatedSamples(100_000, variables, correlation, 'seed-123');

const revenue = samples.map(row => row[0]);
const stats = computeStats(revenue);
console.log(stats.p10, stats.p50, stats.p90);

Latin Hypercube + Iman–Conover (pro mode)

const samples = generateCorrelatedSamples(50_000, variables, correlation, 'seed-123', {
  method: 'lhs',                 // stratified sampling — same accuracy, ~10x fewer sims
  copula: 'iman-conover',        // rank reordering: marginals preserved EXACTLY
  correlationInput: 'spearman',  // treat matrix entries as rank correlations
});

With plain Monte Carlo, 10,000 sims of a normal(100, 10) gives you a mean somewhere in 100 ± 0.2. With LHS the strata force one sample per probability slice, so low-order moments converge dramatically faster — the same trick @RISK and Crystal Ball sell as a headline feature.

Finance statistics

import { expectedShortfall, probabilityAbove, quantile, tailMean } from '@veenie/risk';

expectedShortfall(pnl, 0.05);      // CVaR₉₅: average of the worst 5% of outcomes
probabilityAbove(npv, 0);          // P(NPV > 0)
quantile(pnl, 0.01);               // 1st percentile (VaR₉₉, sign convention yours)
tailMean(returns, 0.10, 'upper');  // average of the best decile

Sensitivity / tornado charts

import { spearmanSensitivity, contributionToVariance } from '@veenie/risk';

// Rank correlation of each input against the output
spearmanSensitivity(inputSamples, outputSamples, [
  { id: 'rev', label: 'Revenue' },
  { id: 'cogs', label: 'COGS' },
]);

// Squared-rho normalized to 100% — the "% of variance explained" view
contributionToVariance(inputSamples, outputSamples, meta);

CLI

The package ships a CLI for quick sanity checks and piping stats into scripts:

# PERT cost estimate, 50k sims
npx @veenie/risk --dist pert --base 100 --min 50 --mode 100 --max 150 --sims 50000

# Latin Hypercube instead of plain Monte Carlo
npx @veenie/risk --dist lognormal --base 100 --spread 0.3 --sims 10000 --method lhs

# Heavy-tailed: Pareto with α=3, reproducible seed
npx @veenie/risk --dist pareto --base 100 --alpha 3 --seed q3-forecast --sims 100000

# Discrete: how many defects in 500 units at 2% failure rate?
npx @veenie/risk --dist binomial --n 500 --p 0.02 --base 500

# Student-t for fat-tailed returns (df=5)
npx @veenie/risk --dist student-t --base 100 --spread 0.15 --df 5

# JSON output for piping
npx @veenie/risk --dist gamma --base 100 --spread 0.5 --json | jq .p90

Flags: --dist --base --sims --seed --method (random|lhs) --json, plus whatever params the distribution takes (--spread --min --max --mode --lambda --p --n --r --df --alpha --shape --mean --center --concentration).

The 17 distributions

Every distribution is base-anchored: you give it a base value (your spreadsheet cell) and shape params, and the distribution centers its mean on the base. This is the core MonteSheet UX contract — swap a normal for a lognormal and your expected value doesn't move.

| Tier | Distribution | Params | Use it for | |------|-------------|--------|-----------| | free | normal | spread (CV) | the default; symmetric uncertainty | | free | triangular | min mode max | quick 3-point estimates | | free | uniform | min max | "anywhere in this range" | | free | discrete-uniform | min max | integer ranges (units, headcount) | | pro | lognormal | spread (CV) | prices, durations — positive & right-skewed | | pro | pert | min mode max | expert 3-point estimates, smoother than triangular | | pro | beta | center concentration | rates & percentages on [0,1] | | pro | gamma | spread (CV) | waiting times, insurance claims | | pro | weibull | shape | failure rates, wind speeds | | pro | exponential | mean | time between events | | pro | student-t | spread df | fat-tailed returns | | pro | pareto | alpha | power laws: wealth, losses, city sizes | | pro | poisson | lambda | event counts | | pro | binomial | n p | successes in n trials | | pro | bernoulli | p | single yes/no events | | pro | geometric | p | trials until first success | | pro | negative-binomial | r p | failures before r-th success |

Each one implements sample, pdf, cdf, and icdf — the CDF being available on everything means P(X > target) queries work uniformly.

What's under the hood

• RNG: Mulberry32 with djb2 string-seed hashing. Same seed → bit-identical results, forever. No global state anywhere (v0.1 had a Box–Muller cache that broke this — see changelog).
• Special functions (src/math/special.ts): Lanczos log-gamma, regularized incomplete gamma (series + continued fraction) and incomplete beta to ~1e-15. Normal CDF via erfc is accurate to machine precision including relative accuracy deep in the tails (Φ(−8) good to 6+ digits), normal ICDF is Acklam's algorithm plus one Halley refinement step.
• Correlation: Gaussian copula (NORTA) via Cholesky, with automatic repair of non-positive-definite matrices (shrink toward identity, warn with the shrinkage %). Or Iman–Conover rank reordering, which permutes rather than transforms — marginals and LHS stratification survive exactly.
• Spearman↔Pearson: correlationInput: 'spearman' converts via the exact Gaussian-copula relation ρ = 2·sin(πρₛ/6).
• Tests: 140 tests, 40k+ assertions. Every distribution is pinned to closed-form analytic moments, round-tripped through CDF∘ICDF on a 999-point grid, and its PDF numerically integrated back to 1.

Quirks & conventions (read before shipping on top of this)

• stdDev is population (÷n), matching computeStats since 0.1. Use sampleStdDev for Excel STDEV.S (÷n−1) semantics.
• Discrete quantiles use the standard right-continuous inverse: icdf(p) returns the smallest integer k with CDF(k) ≥ p.
• Probabilities are clamped to [1e-12, 1−1e-12] before inverse transforms, so unbounded distributions never return ±Infinity from the sampler. normalICDF(0) directly still returns −Infinity, as it should.
• Seed streams changed in 0.2.0. Same seed gives different (but distributionally identical) samples than 0.1.x, because the copula path now consumes uniforms via inverse transform instead of polar sampling. Statistical results are equivalent; exact sample values are not. Pin 0.1.x if you froze golden outputs.
• Non-PD correlation matrices are repaired, not dropped. 0.1.x silently fell back to independent sampling; 0.2.0 shrinks the matrix toward identity just enough to factor, and warns.
• Correlation matrices are symmetrized (upper/lower averaged) and off-diagonals clamped to ±0.9999 before use.
• Base-anchoring edge cases: pareto requires α > 1 (else the mean doesn't exist — the param floor enforces this), student-t requires df > 2 for the spread to mean what you think it means.

Why this exists

If you build financial models, project cost estimators, or engineering tolerance sims in JS/TS, the alternative is stitching together five unmaintained math libraries with inconsistent parameterizations. @veenie/risk puts the hard parts — inverse CDF sampling, matrix decomposition, correlation induction, tail-accurate special functions — behind one tested, seeded, tree-shakeable API. It feeds directly into headless calc engines like HyperFormula, web workers, or interactive calculators.

Want it inside your spreadsheet instead of your codebase? That's MonteSheet — =MS_NORMAL(100, 0.1) in a cell and go.

License

MIT. See LICENSE. Built by Veenie.