Monte Carlo Cost-to-Complete Simulation in Python

A single-number forecast-to-complete answers the wrong question for a completion guarantor. “The remaining shoot will cost 2.4 million” tells an underwriter nothing about how badly it could go, and the whole point of a contingency is to absorb the bad days that a point estimate pretends away. What a bond covenant actually needs is a distribution: given the variability in daily shoot cost, how much cost-to-complete would cover the project 80 percent of the time, or 95 percent of the time? A Monte Carlo simulation answers exactly that by sampling many possible completions and reading the percentiles off the result. The catch is that a simulation whose number changes every time you run it is worthless to an auditor, so this walkthrough builds a seeded, reproducible simulation over remaining shoot days, keeps money as Decimal at the reporting boundary, and emits P50, P80, and P95 bounds a guarantor can recompute to the same figures months later. It extends the point-estimate arithmetic in Cost-to-Complete Projections with a bounded, probabilistic view of the soft forecast.

Prerequisites and Context

This page assumes a validated cost-to-complete baseline already exists — the actuals, open commitments, and point forecast from the parent Cost-to-Complete Projections engine — and focuses on turning the remaining uncommitted work into a distribution. It targets Python 3.11+ and a deliberately small stack: the standard-library random module seeded explicitly for reproducibility, statistics for percentile extraction, decimal for money at the reporting boundary, hashlib for the run fingerprint, and zoneinfo for the timezone-aware timestamp; Pydantic v2 (model_validate, field_validator) validates the simulation parameters. The simulation samples in float for speed inside the loop and crosses back to Decimal only at the boundary where results become reported money — the one place the conversion is allowed, and only through Decimal(str(x)). The remaining-day parameters — expected daily cost, its variability, and the number of shoot days left — are derived from the same cost-coded baseline standardized through Cost Code Standardization, and the bounded output feeds the reserve logic in Contingency Drawdown Tracking.

Step-by-Step Implementation

The simulation proceeds in four steps. First, validate the parameters and pin an explicit integer seed. Second, run many trials, each summing a sampled per-day cost across the remaining shoot days to produce one possible cost-to-complete. Third, read P50, P80, and P95 off the sorted trials. Fourth, cross from float to Decimal at that boundary, size the contingency as the gap between the point forecast and the chosen percentile, and emit a hashed audit record that pins the seed. The seed is the whole ballgame: a fixed seed makes random produce the identical sequence of draws on every run, which is what converts a stochastic method into a deterministic, auditable one.

from __future__ import annotations

import hashlib
import json
import random
from datetime import datetime
from decimal import Decimal, ROUND_HALF_UP
from statistics import quantiles
from zoneinfo import ZoneInfo

from pydantic import BaseModel, ConfigDict, Field, field_validator

CENTS = Decimal("0.01")


class SimParams(BaseModel):
    """Parameters for a seeded cost-to-complete simulation."""

    model_config = ConfigDict(frozen=True)

    remaining_days: int = Field(gt=0)
    mean_daily_cost: Decimal = Field(gt=Decimal("0"))
    daily_cost_stddev: Decimal = Field(ge=Decimal("0"))
    point_forecast: Decimal = Field(ge=Decimal("0"))  # deterministic ETC baseline
    trials: int = Field(gt=0, le=1_000_000)
    seed: int  # pinned for reproducibility — same seed, same bounds

    @field_validator(
        "mean_daily_cost", "daily_cost_stddev", "point_forecast",
        mode="before",
    )
    @classmethod
    def reject_float(cls, v) -> Decimal:
        if isinstance(v, float):
            raise ValueError("monetary fields must be str/Decimal, never float")
        return Decimal(str(v))


def _q(amount: Decimal) -> Decimal:
    return amount.quantize(CENTS, rounding=ROUND_HALF_UP)


def simulate_cost_to_complete(p: SimParams) -> dict:
    # Seed the module RNG explicitly. A pinned seed => an identical draw
    # sequence => identical percentiles on every re-run.
    rng = random.Random(p.seed)
    mean = float(p.mean_daily_cost)
    sigma = float(p.daily_cost_stddev)

    outcomes: list[float] = []
    for _ in range(p.trials):
        total = 0.0
        for _ in range(p.remaining_days):
            # Truncate at zero: a shoot day cannot cost negative money.
            day = rng.gauss(mean, sigma)
            total += day if day > 0.0 else 0.0
        outcomes.append(total)

    outcomes.sort()
    # 100 cut points => percentile buckets; index i is the (i+1)th percentile.
    cuts = quantiles(outcomes, n=100, method="inclusive")
    p50 = _q(Decimal(str(cuts[49])))
    p80 = _q(Decimal(str(cuts[79])))
    p95 = _q(Decimal(str(cuts[94])))

    # Contingency sized to cover the P80 outcome above the point forecast.
    contingency_p80 = _q(max(p80 - p.point_forecast, Decimal("0")))
    contingency_p95 = _q(max(p95 - p.point_forecast, Decimal("0")))

    payload = {
        "remaining_days": p.remaining_days,
        "trials": p.trials,
        "seed": p.seed,
        "mean_daily_cost": str(_q(p.mean_daily_cost)),
        "daily_cost_stddev": str(_q(p.daily_cost_stddev)),
        "point_forecast": str(_q(p.point_forecast)),
        "p50_cost_to_complete": str(p50),
        "p80_cost_to_complete": str(p80),
        "p95_cost_to_complete": str(p95),
        "contingency_to_p80": str(contingency_p80),
        "contingency_to_p95": str(contingency_p95),
    }
    canonical = json.dumps(payload, sort_keys=True, separators=(",", ":")).encode()
    payload["payload_sha256"] = hashlib.sha256(canonical).hexdigest()
    payload["computed_at"] = datetime.now(tz=ZoneInfo("UTC")).isoformat()
    return payload

Because the run is driven by an explicitly seeded random.Random(p.seed) instance rather than the global, shared generator, it is isolated and reproducible: the same parameters and the same seed produce the same sorted outcomes, the same percentiles, and the same hash on every machine. Two properties make the number defensible. The seed is a first-class, logged parameter — not an implementation detail — so an auditor reproducing the run supplies the same integer and gets the same bounds. And money crosses from float to Decimal only at the percentile boundary, through Decimal(str(...)), so the fast inner loop never contaminates a reported figure with a value that was already Decimal upstream.

Seeded Monte Carlo cost-to-complete simulation flow Validated simulation parameters and a pinned integer seed enter a seeded random generator. The generator runs many trials, each summing a sampled per-day cost across the remaining shoot days into one possible cost-to-complete outcome. The sorted outcomes are read at three percentiles — P50, P80, and P95 — which cross from float to Decimal at the reporting boundary. The chosen percentile above the point forecast sizes the contingency, and the whole run is serialized with its seed into a SHA-256 stamped audit record so it reproduces exactly. Params + pinned seed days · mean · stddev Seeded trials sum sampled day cost × remaining days Sort → P50 · P80 · P95 float → Decimal boundary Size contingency percentile − point forecast Serialize run + seed + SHA-256 audit hash same seed → same bounds

Audit Trail Requirements

A simulation is only an audit instrument if a third party can reproduce it exactly, and reproduction requires that the record carry everything the run consumed. At minimum the persisted payload must log: the seed, the trials count, and the remaining_days; the mean_daily_cost and daily_cost_stddev that shaped the distribution; the point_forecast baseline; the three reported bounds p50_cost_to_complete, p80_cost_to_complete, and p95_cost_to_complete as Decimal strings; the contingency_to_p80 and contingency_to_p95 sizes; and the payload_sha256 fingerprint over the canonical payload. The seed and the trial count are not optional metadata — they are the reproduction key, and a bound reported without them is a number no guarantor can verify. Persist the record to append-only, write-once storage before any contingency decision commits against it, and treat a re-simulation with revised parameters as a new report-dated record rather than an edit, so the bound a covenant test was measured against stays recoverable. This mirrors the snapshot discipline of the parent Cost-to-Complete Projections engine, where a report-dated, hashed artifact is the unit of truth.

Gotchas and Production Edge Cases

Seed pinning and generator isolation. The single failure that voids the whole exercise is an unpinned or shared generator. Calling the module-level random.gauss without seeding, or seeding the global generator that other code also draws from, makes the run non-reproducible — a guarantor who re-runs it gets different bounds and rightly distrusts the report. Always construct a dedicated random.Random(seed) instance, log the seed, and never let unrelated code share it. If you later move to a numeric library’s generator for speed, pin its seed the same way and record which generator produced the run, because different generators yield different sequences from the same seed.

Float-to-Decimal boundary discipline. The inner loop runs in float deliberately — sampling and summing millions of draws in Decimal is needlessly slow — but every value that becomes reported money must cross to Decimal through its string form at the percentile boundary, exactly once. Never build a Decimal from a raw float with Decimal(x); that imports the binary rounding error the boundary exists to stop. And never re-multiply a point_forecast that arrived as Decimal from the upstream projection — it is already exact, so mixing it into the float loop would corrupt it. Keep the float world inside the loop and the Decimal world at the edges, with Decimal(str(...)) as the only bridge.

Choosing the covenant percentile. P50 is a coin flip and no basis for a reserve; a bond covenant that wants genuine protection sizes contingency to P80 or P95, and which one is a negotiated term, not an engineering default. Report all three so the guarantor can choose, and size the contingency as the gap between the chosen percentile and the point forecast rather than the percentile itself — the point forecast is already funded in the estimate-at-completion, and only the tail above it needs reserve. A contingency sized to the raw P95 double-counts the money already in the baseline.

Distribution realism. A symmetric normal draw understates real production risk, where the downside tail is fatter than the upside — a lost location or a weather day costs far more than a good day saves. Truncating at zero, as the loop does, prevents the physical absurdity of a negative-cost day, but a production-grade model often reaches for a right-skewed distribution to reflect that overruns cluster on the expensive side. Whatever the shape, it must be a logged parameter of the run so the audit record captures which distribution produced the bound.

Idempotency on replay. Because the run is seeded and the store is append-only, re-running after a partial failure is safe only if the consumer deduplicates on (seed, trials, remaining_days, payload_sha256). A retried batch with the identical seed reproduces the identical bounds, so the guard collapses the duplicate rather than booking a second contingency figure for the same simulation. When an upstream baseline is missing or delayed, route the parameter set through the controlled path in Contingency Drawdown Tracking rather than simulating against a guessed forecast inline.

Up one level: Cost-to-Complete Projections, part of Completion Bond Reporting & Guarantor Analytics.