Supply Chain Risk & Interdependent Security

Your security posture is not determined by your controls alone. It is determined by the intersection of your controls and your vendors’ controls, your vendors’ vendors’ controls, and the shared infrastructure underneath all of it.

This is the problem of interdependent security (Kunreuther & Heal, 2003): each organization’s risk depends on the security investments of others in the network. You can have a perfect patch cadence, world-class SOC, and zero known vulnerabilities – and still get breached because your HVAC vendor reuses passwords.

Supply chain risk is not additive. It is not even linear. It exhibits correlation, contagion, and cascade dynamics that make naive risk calculations dangerously wrong.

This notebook builds up from the simplest model (independent vendor failures) to correlated failures, cascade propagation, and investment.

Section 6 repeats the exercise on a real software dependency graph.

Setup

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from decision_security.synth import make_rng, sample
from decision_security.montecarlo import simulate_aggregate_losses, make_lognormal_severity, var_es
from decision_security.survival import km_estimator
from decision_security.viz import plot_km

rng = make_rng(42)

PRIMARY = "#1A1A1A"
ACCENT = "#E74C3C"
DARK_BG = "#34495E"
LIGHT_GRAY = "#95A5A6"

plt.rcParams.update({
    "font.family": "serif",
    "font.size": 10,
    "axes.labelsize": 11,
    "axes.titlesize": 12,
    "xtick.labelsize": 9,
    "ytick.labelsize": 9,
    "legend.fontsize": 9,
    "figure.dpi": 150,
    "axes.spines.top": False,
    "axes.spines.right": False,
    "text.color": PRIMARY,
    "axes.labelcolor": PRIMARY,
    "axes.edgecolor": PRIMARY,
    "xtick.color": PRIMARY,
    "ytick.color": PRIMARY,
})

1. Your Risk Is Not Just Yours

Organization A depends on three vendors: B, C, and D. Each vendor has an independent probability of suffering a breach that exposes A’s data or systems. A breach at any vendor compromises A.

Under independence, A’s compound exposure is one minus the probability that every vendor survives:

\[P(\text{A exposed}) = 1 - \prod_{i} (1 - p_i)\]

This is the complement of the joint survival probability under independence – each factor \((1-p_i)\) is a vendor surviving.

The growth is sublinear (concave) but visible. Each vendor you add multiplies the survival probability by another factor below 1, so exposure keeps climbing toward 1 – but in shrinking increments. Adding a fourth vendor at “only” 3% still ratchets A’s exposure up by 2.3 points – less than 3 points, because each new vendor acts only on the shrinking probability that everyone else survived. The danger is not that any one step is large; it is that the steps never stop and never reverse.

# --- Compound vendor exposure under independence ---

vendor_probs = {"Vendor B": 0.05, "Vendor C": 0.08, "Vendor D": 0.12}
probs = list(vendor_probs.values())
compound_3 = 1.0 - np.prod([1 - p for p in probs])

print("Individual vendor breach probabilities:")
for name, p in vendor_probs.items():
    print(f"  {name}: {p:.0%}")
print(f"\nCompound exposure (3 vendors): {compound_3:.1%}")
print(f"Sum of individual probs:        {sum(probs):.1%}")
print("Below the sum (overlap), but far above any single vendor's risk.")

vendor_probs_4 = {**vendor_probs, "Vendor E": 0.03}
probs_4 = list(vendor_probs_4.values())
compound_4 = 1.0 - np.prod([1 - p for p in probs_4])

print(f"\n--- Adding Vendor E (3% breach probability) ---")
print(f"Compound exposure (4 vendors): {compound_4:.1%}")
print(f"Marginal increase from Vendor E: {compound_4 - compound_3:.1%}")

Individual vendor breach probabilities:
  Vendor B: 5%
  Vendor C: 8%
  Vendor D: 12%

Compound exposure (3 vendors): 23.1%
Sum of individual probs:        25.0%
Below the sum (overlap), but far above any single vendor's risk.

--- Adding Vendor E (3% breach probability) ---
Compound exposure (4 vendors): 25.4%
Marginal increase from Vendor E: 2.3%

# --- Bar chart: individual vendor risk vs compound exposure ---

fig, axes = plt.subplots(1, 2, figsize=(10, 4.5))

ax = axes[0]
labels_3 = list(vendor_probs.keys()) + ["A's Exposure"]
values_3 = list(vendor_probs.values()) + [compound_3]
colors_3 = [DARK_BG] * 3 + [ACCENT]
bars = ax.bar(labels_3, values_3, color=colors_3, edgecolor="white")
for bar, v in zip(bars, values_3):
    ax.text(bar.get_x() + bar.get_width() / 2, v + 0.005,
            f"{v:.1%}", ha="center", va="bottom", fontsize=9, fontweight="bold")
ax.set_ylabel("P(breach)")
ax.set_title("3 Vendors")
ax.set_ylim(0, 0.35)

ax = axes[1]
labels_4 = list(vendor_probs_4.keys()) + ["A's Exposure"]
values_4 = list(vendor_probs_4.values()) + [compound_4]
colors_4 = [DARK_BG] * 4 + [ACCENT]
bars = ax.bar(labels_4, values_4, color=colors_4, edgecolor="white")
for bar, v in zip(bars, values_4):
    ax.text(bar.get_x() + bar.get_width() / 2, v + 0.005,
            f"{v:.1%}", ha="center", va="bottom", fontsize=9, fontweight="bold")
ax.set_ylabel("P(breach)")
ax.set_title("4 Vendors (adding E at 3%)")
ax.set_ylim(0, 0.35)

fig.suptitle("Supply chain risk accumulates sublinearly but relentlessly in vendor count",
             fontsize=13, fontweight="bold", y=1.02)
fig.tight_layout()
plt.show()

print("\n--- Scaling: N vendors each with 5% breach probability ---")
print(f"{'N':>3}  {'P(A exposed)':>14}  {'marginal +':>12}")
prev = 0.0
for n in range(1, 21):
    compound = 1 - (1 - 0.05) ** n
    marginal = compound - prev
    if n in (1, 2, 3, 5, 10, 15, 20):
        print(f"{n:>3}  {compound:>14.2%}  {marginal:>12.2%}")
    prev = compound
print("\nThe marginal shrinks (5.00% -> 4.75% -> 4.51% ...) but never turns negative.")

--- Scaling: N vendors each with 5% breach probability ---
  N    P(A exposed)    marginal +
  1           5.00%         5.00%
  2           9.75%         4.75%
  3          14.26%         4.51%
  5          22.62%         4.07%
 10          40.13%         3.15%
 15          53.67%         2.44%
 20          64.15%         1.89%

The marginal shrinks (5.00% -> 4.75% -> 4.51% ...) but never turns negative.

2. Correlated Failures

The independence assumption in Section 1 is convenient but wrong. Vendors share infrastructure (AWS, Azure), run the same software (Log4j, MOVEit, SolarWinds Orion have shown this effect), and face the same threat actors; their risk also gets misjudged by the same rating services.

We model this with a common-cause factor: with probability \(p_{\text{systemic}} = 0.15\), a systemic event occurs that breaches ALL vendors simultaneously. Otherwise, vendors fail independently as before.

Two points about this model. First, it is not mean-preserving: bolting a 15% “everything breaches” event on top of the independent baseline raises both the mean and the tail. Second, that mean jump is an artifact of how we injected the correlation, not of correlation itself. A mean-preserving correlation model (same marginal breach probabilities, just made dependent) would leave the mean unchanged and still fatten the tail. The structural point is the tail: dependence concentrates probability mass on the all-at-once outcome, which is exactly what VaR and ES are built to expose.

# --- Independent vs Correlated vendor failure model ---

N_SIM = 10_000
vendor_p = np.array([0.05, 0.08, 0.12])
vendor_loss = np.array([500_000, 800_000, 1_200_000])  # loss if vendor breached
p_systemic = 0.15

breach_ind = rng.random((N_SIM, len(vendor_p))) < vendor_p
losses_independent = breach_ind @ vendor_loss

# When the systemic event fires, all vendors breach together. That is an extreme
# dependence structure; real supply chains sit between this and independence.
breach_cor = rng.random((N_SIM, len(vendor_p))) < vendor_p
losses_correlated = breach_cor @ vendor_loss
systemic = rng.random(N_SIM) < p_systemic
losses_correlated[systemic] = vendor_loss.sum()

var95_ind, es95_ind = var_es(losses_independent, alpha=0.95)
var95_cor, es95_cor = var_es(losses_correlated, alpha=0.95)

print(f"{'Metric':<25} {'Independent':>14} {'Correlated':>14}")
print(f"{'':─<25} {'':─>14} {'':─>14}")
print(f"{'Mean loss':<25} ${losses_independent.mean():>12,.0f} ${losses_correlated.mean():>12,.0f}")
print(f"{'Median loss':<25} ${np.median(losses_independent):>12,.0f} ${np.median(losses_correlated):>12,.0f}")
print(f"{'VaR (95%)':<25} ${var95_ind:>12,.0f} ${var95_cor:>12,.0f}")
print(f"{'ES (95%)':<25} ${es95_ind:>12,.0f} ${es95_cor:>12,.0f}")
print(f"{'P(loss > $1M)':<25} {(losses_independent > 1_000_000).mean():>14.1%} {(losses_correlated > 1_000_000).mean():>14.1%}")
print(f"{'P(total wipeout)':<25} {(losses_independent >= vendor_loss.sum()).mean():>14.2%} {(losses_correlated >= vendor_loss.sum()).mean():>14.2%}")
print(f"\nMean loss rises {losses_correlated.mean() / losses_independent.mean():.1f}x "
      f"(the systemic event is not mean-preserving),")
print(f"VaR95 rises {var95_cor / max(var95_ind, 1):.1f}x, and ES95 rises "
      f"{es95_cor / max(es95_ind, 1):.1f}x.")

Metric                       Independent     Correlated
───────────────────────── ────────────── ──────────────
Mean loss                 $     230,750 $     570,780
Median loss               $           0 $           0
VaR (95%)                 $   1,200,000 $   2,500,000
ES (95%)                  $   1,436,727 $   2,500,000
P(loss > $1M)                      12.2%          25.4%
P(total wipeout)                   0.04%         14.98%

Mean loss rises 2.5x (the systemic event is not mean-preserving),
VaR95 rises 2.1x, and ES95 rises 1.7x.

# --- Overlaid histograms: independent vs correlated ---
from matplotlib.ticker import FuncFormatter

usd_k = FuncFormatter(lambda x, _: f"${x / 1e3:,.0f}K")

fig, ax = plt.subplots(figsize=(9, 5))

bins = np.linspace(0, vendor_loss.sum() * 1.05, 50)

# Filled bars for one series and an outline for the other, so no blended colour appears.
ax.hist(losses_correlated, bins=bins, color=ACCENT,
        label="Common-cause model (15% systemic)", edgecolor="white", linewidth=0.3)
ax.hist(losses_independent, bins=bins, histtype="step", color=DARK_BG, linewidth=1.5,
        label="Independent model")

ax.axvline(var95_ind, color=DARK_BG, linestyle="--", linewidth=1.5)
ax.axvline(var95_cor, color=ACCENT, linestyle="--", linewidth=1.5)

ymax = ax.get_ylim()[1]
ax.text(var95_ind - 30_000, ymax * 0.5, f"VaR95\n${var95_ind:,.0f}",
        fontsize=8, color=DARK_BG, ha="right", va="bottom")
ax.text(var95_cor - 30_000, ymax * 0.5, f"VaR95\n${var95_cor:,.0f}",
        fontsize=8, color=ACCENT, ha="right", va="bottom")

ax.xaxis.set_major_formatter(usd_k)
ax.set_xlabel("Total supply chain loss")
ax.set_ylabel("Frequency")
ax.set_title("Supply chain loss distribution: independent vs common-cause model (15% systemic)")
ax.legend(loc="upper right", frameon=True, fancybox=False, edgecolor=LIGHT_GRAY)

fig.tight_layout()
plt.show()

print("The spike at $2.5M is the systemic event: all three vendors breached at once, in ~15%")
print("of runs, and all but absent from the independent model. It lifts the mean AND stretches")
print("the tail; a mean-preserving dependence model would move only the tail, the part that")
print("changes capital and insurance decisions.")

The spike at $2.5M is the systemic event: all three vendors breached at once, in ~15%
of runs, and all but absent from the independent model. It lifts the mean AND stretches
the tail; a mean-preserving dependence model would move only the tail, the part that
changes capital and insurance decisions.

Check it in your own environment.

1. List what your critical vendors share: cloud provider and region, key software components, identity provider, managed service partners.

2. For each shared item, ask which vendors would fail together if it failed, and model that as one scenario rather than separate ones.

3. Check whether your aggregate loss figures add vendor risks as if they were independent.

3. Cascade Dynamics

Beyond correlation, breaches can propagate: a compromised vendor becomes a launch point for attacking its customers and partners. SolarWinds was not just a correlated failure: the trojanized Orion update reached about 18,000 customers, and at Mimecast the attackers used that foothold to steal a certificate and reach Mimecast’s own customers.

We model 10-node networks where each node is connected to several others. A breach at any node propagates to each connected neighbor with probability \(p_{\text{propagate}} = 0.3\) per edge. Starting from one randomly chosen compromised node, we iterate until no new nodes are compromised. To avoid drawing conclusions from one lucky graph, we average over 20 random graphs with random seed nodes.

Because the network has only 10 nodes, the cascade size is bounded – so the distribution is not heavy-tailed in the extreme-value sense. What it is is right-skewed and high-variance: most cascades fizzle at 1-2 nodes (the seed’s edges happen to fail), but a substantial minority escalate, and a few sweep most of the network. The mean (2.5) is close to the median (2) and says nothing about the 15% of cascades that reach 5+ nodes or the 2.6% that reach 8+.

# --- Cascade simulation over many random 10-node networks ---

N_NODES = 10
P_PROPAGATE = 0.3
EDGE_DENSITY = 0.3
N_GRAPHS = 20
SIMS_PER_GRAPH = 50

graph_rng = make_rng(99)      # separate stream, so graph structure does not shift the cascade draws


def random_graph(n_nodes, density, rng):
    """Erdos-Renyi-style undirected graph as an adjacency list."""
    adj = {i: [] for i in range(n_nodes)}
    for i in range(n_nodes):
        for j in range(i + 1, n_nodes):
            if rng.random() < density:
                adj[i].append(j)
                adj[j].append(i)
    return adj


def simulate_cascade(adj, start_node, p_propagate, rng):
    """BFS cascade: each compromised node infects each neighbor with probability p."""
    compromised = {start_node}
    frontier = [start_node]
    while frontier:
        next_frontier = []
        for node in frontier:
            for neighbor in adj[node]:
                if neighbor not in compromised and rng.random() < p_propagate:
                    compromised.add(neighbor)
                    next_frontier.append(neighbor)
        frontier = next_frontier
    return len(compromised)


cascade_sizes = []
edge_counts = []
for _ in range(N_GRAPHS):
    adj = random_graph(N_NODES, EDGE_DENSITY, graph_rng)
    edge_counts.append(sum(len(v) for v in adj.values()) // 2)
    for _ in range(SIMS_PER_GRAPH):
        start = int(rng.integers(0, N_NODES))
        cascade_sizes.append(simulate_cascade(adj, start, P_PROPAGATE, rng))
cascade_sizes = np.array(cascade_sizes)

print(f"Averaged over {N_GRAPHS} random graphs ({N_NODES} nodes, "
      f"~{np.mean(edge_counts):.0f} edges each), {SIMS_PER_GRAPH} random-seed cascades per graph")
print(f"\nCascade size distribution ({len(cascade_sizes):,} simulations):")
print(f"  Mean:   {cascade_sizes.mean():.1f} nodes")
print(f"  Median: {np.median(cascade_sizes):.0f} nodes")
print(f"  P(== 1 node, fizzles): {(cascade_sizes == 1).mean():.1%}")
print(f"  P(>= 5 nodes):         {(cascade_sizes >= 5).mean():.1%}")
print(f"  P(>= 8 nodes):         {(cascade_sizes >= 8).mean():.1%}")
print(f"  Max cascade:           {cascade_sizes.max()} nodes")

Averaged over 20 random graphs (10 nodes, ~12 edges each), 50 random-seed cascades per graph

Cascade size distribution (1,000 simulations):
  Mean:   2.5 nodes
  Median: 2 nodes
  P(== 1 node, fizzles): 48.0%
  P(>= 5 nodes):         15.3%
  P(>= 8 nodes):         2.6%
  Max cascade:           10 nodes

# --- Distribution of cascade sizes ---

fig, ax = plt.subplots(figsize=(8, 5))

counts, edges, bars = ax.hist(cascade_sizes, bins=np.arange(0.5, N_NODES + 1.5, 1),
                               color=DARK_BG, edgecolor="white", linewidth=0.5)

for i, bar in enumerate(bars):
    if edges[i] + 0.5 >= 5:  # 5+ nodes = tail event
        bar.set_color(ACCENT)

ax.axvline(cascade_sizes.mean(), color=ACCENT, linestyle="--", linewidth=1.5,
           label=f"Mean = {cascade_sizes.mean():.1f}")
ax.axvline(np.median(cascade_sizes), color=LIGHT_GRAY, linestyle=":", linewidth=1.5,
           label=f"Median = {np.median(cascade_sizes):.0f}")

ax.set_xlabel("Cascade size (nodes compromised)")
ax.set_ylabel("Frequency")
ax.set_title("Cascade size distribution (20 graphs x 50 seeds = 1,000 simulations)")
ax.set_xticks(range(1, N_NODES + 1))
ax.legend(loc="upper right", frameon=True, fancybox=False, edgecolor=LIGHT_GRAY)

tail_pct = (cascade_sizes >= 5).mean()
ax.annotate(f"{tail_pct:.0%} of cascades\nreach 5+ nodes",
            xy=(5, counts[4] if len(counts) > 4 else 0),
            xytext=(7, max(counts) * 0.6),
            fontsize=9, color=ACCENT,
            arrowprops=dict(arrowstyle="->", color=ACCENT, lw=1.2))

fig.tight_layout()
plt.show()

print("Right-skewed and high-variance: most cascades fizzle near the seed, a shoulder spreads")
print("partway, and a thin tail takes most of the network. The mean "
      f"({cascade_sizes.mean():.1f}) is close to the median ({np.median(cascade_sizes):.0f}),")
print("so the average case hides the tail. The upper tail is what justifies the investment.")

Right-skewed and high-variance: most cascades fizzle near the seed, a shoulder spreads
partway, and a thin tail takes most of the network. The mean (2.5) is close to the median (2),
so the average case hides the tail. The upper tail is what justifies the investment.

4. Can You Tell Whether Vendor Assessments Work?

Vendor assessments are the standard response to supply chain risk, but their evidence value is weak. Most are questionnaires with self-reported yes/no answers. The largest providers do not answer custom questionnaires; they point to pre-completed standard ones (CAIQ, SIG), audit reports such as SOC 2 (an independent attestation), and a shared-responsibility model. Audit reports attest that controls exist and operated; they say little about how likely a breach is.

So does an assessment program delay a vendor-originated compromise? We do not know, and the question here is whether you could know. Suppose a program really did cut the hazard of compromise by some amount. If you compared organizations with and without the program over a 36-month window, how often would the difference show up in the data?

We simulate that comparison many times, for different effect sizes and numbers of organizations, and count how often a standard log-rank test detects the effect.

# --- If a program had a real effect, how often would the data show it? ---
from scipy.stats import logrank, CensoredData

CENSOR_AT = 36        # months -- observation window
# Assumed baseline without a program: Weibull k=1.2 (hazard rises with time), lambda=14,
# so about 96% of organizations have a vendor-originated compromise within 36 months.
# That is a best case for detection; power falls sharply at lower, more realistic rates.
K, LAM = 1.2, 14.0
N_REPS = 200          # simulated studies per cell

def simulate_group(n, hazard_ratio, rng):
    # Scaling a Weibull hazard by HR stretches its scale by HR**(-1/k).
    t = sample("weibull", n, rng=rng, k=K, lam=LAM * hazard_ratio ** (-1 / K))
    return np.minimum(t, CENSOR_AT), t <= CENSOR_AT

def detection_rate(n, hazard_ratio, rng):
    hits = 0
    for _ in range(N_REPS):
        t0, e0 = simulate_group(n, 1.0, rng)
        t1, e1 = simulate_group(n, hazard_ratio, rng)
        p = logrank(CensoredData.right_censored(t0, ~e0), CensoredData.right_censored(t1, ~e1)).pvalue
        hits += p < 0.05
    return hits / N_REPS

power_rng = make_rng(7)
effects = {"hazard cut 10%": 0.9, "cut 20%": 0.8, "cut 50%": 0.5}
sizes = [50, 200, 500]
power = pd.DataFrame({label: [detection_rate(n, hr, power_rng) for n in sizes]
                      for label, hr in effects.items()},
                     index=[f"{n} orgs per group" for n in sizes])
print(f"Assumed baseline: {1 - np.exp(-(CENSOR_AT / LAM) ** K):.0%} compromised within {CENSOR_AT} months without a program.")
print("Share of simulated studies in which the log-rank test detects the effect (p < 0.05):\n")
print(power.map(lambda v: f"{v:.0%}").to_string())

Assumed baseline: 96% compromised within 36 months without a program.
Share of simulated studies in which the log-rank test detects the effect (p < 0.05):

                   hazard cut 10% cut 20% cut 50%
50 orgs per group              8%     24%     92%
200 orgs per group            20%     64%    100%
500 orgs per group            34%     93%    100%

# --- One 50-vs-50 comparison with a real 20% effect ---
one_rng = make_rng(3)
t0, e0 = simulate_group(50, 1.0, one_rng)
t1, e1 = simulate_group(50, 0.8, one_rng)
p_one = logrank(CensoredData.right_censored(t0, ~e0), CensoredData.right_censored(t1, ~e1)).pvalue

fig, ax = plt.subplots(figsize=(8, 5))
for (t, e), label, col in [((t0, e0), "Without program", LIGHT_GRAY),
                           ((t1, e1), "With program (true hazard cut 20%)", DARK_BG)]:
    km_t, km_s = km_estimator(t, e.astype(int))
    km_t, km_s = np.r_[0, km_t], np.r_[1.0, km_s]  # start the curve at (0, 1)
    ax.step(km_t, km_s, where="post", color=col, linewidth=2, label=label)
ax.set_xlabel("Time (months)")
ax.set_ylabel("Proportion not yet compromised")
ax.set_title(f"One study, 50 orgs per group: log-rank p = {p_one:.2f}")
ax.set_xlim(0, CENSOR_AT + 1)
ax.set_ylim(0, 1.05)
ax.legend(loc="upper right", frameon=False)
fig.tight_layout()
plt.show()

print(f"A real 20% effect is detected in {power.loc['50 orgs per group', 'cut 20%']:.0%} of 50-vs-50 studies "
      f"and {power.loc['200 orgs per group', 'cut 20%']:.0%} of 200-vs-200 studies.")

A real 20% effect is detected in 24% of 50-vs-50 studies and 64% of 200-vs-200 studies.

Unless the effect is large, a single organization cannot tell from breach outcomes whether its assessment program works, and a 20% effect is missed in about a third of 200-vs-200 studies even under this generous base rate. The table is a best case: it assumes about 96% of organizations without a program are compromised within 36 months, and power falls sharply at more realistic, lower base rates. And a real with-vs-without comparison is observational: a difference, if found, would still not show the program caused it. Claims that assessments help, and claims that they are useless, are both ahead of the evidence. Controls on your own side of the vendor connection (least-privilege vendor access, segmented connections, monitoring of vendor-originated activity) at least let you verify that they are in place; their effect on breach outcomes is still an estimate.

5. Comparing Vendor Security Investments

You have $100K to invest in vendor security. Three strategies:

Strategy	Description	Mechanism
(a) Spread audit	$25K per vendor (4 vendors)	Reduces each vendor’s P(breach) by 20%
(b) Deep audit	$100K on highest-risk vendor (D)	Reduces vendor D’s P(breach) by 60%
(c) Contractual controls	$100K on SLAs and contracts	Assumes contractual recovery offsets 50% of each breach loss

The reductions in (a) and (b) and the 50% recovery in (c) are assumptions for the exercise – as Section 4 shows, nobody can currently measure how much an audit lowers a vendor’s breach probability. Treat them as scenarios to compare, not expected results.

Strategy (a) is the diversified approach. Strategy (b) is targeted – fix the weakest link. Strategy (c) is transfer – accept breaches will happen, limit the damage.

# --- Three strategies, all costing the same $100K ---

N_SIM = 10_000

vendor_names = ["Vendor B", "Vendor C", "Vendor D", "Vendor E"]
base_probs = np.array([0.05, 0.08, 0.12, 0.03])
base_losses = np.array([500_000, 800_000, 1_200_000, 300_000])

# Common random numbers: one uniform matrix reused for every strategy, so differences
# reflect the strategy, not luck of the draw. Lowering a probability can only remove
# breaches from the identical scenario set.
U = rng.random((N_SIM, len(base_probs)))


def crn_losses(probs, losses, U, recovery=0.0):
    breached = U < probs
    return (breached * (losses * (1 - recovery))).sum(axis=1)


losses_baseline = crn_losses(base_probs, base_losses, U)

# (a) spread audit: assumed 20% reduction in each vendor's P(breach)
losses_a = crn_losses(base_probs * 0.80, base_losses, U)

# (b) deep audit: assumed 60% reduction in vendor D's P(breach)
probs_b = base_probs.copy()
probs_b[2] *= 0.40                # vendor D
losses_b = crn_losses(probs_b, base_losses, U)

# (c) contractual controls: assumed 50% recovery on every vendor breach
losses_c = crn_losses(base_probs, base_losses, U, recovery=0.50)

strategies = {
    "Baseline (no investment)": losses_baseline,
    "(a) Spread audit ($25K x 4)": losses_a,
    "(b) Deep audit ($100K on D)": losses_b,
    "(c) Contracts ($100K)": losses_c,
}

print("All three strategies spend the SAME $100K; figures below are residual breach loss.")
print(f"{'Strategy':<30} {'E[Loss]':>12} {'VaR95':>12} {'ES95':>12}")
print(f"{'':─<30} {'':─>12} {'':─>12} {'':─>12}")
for name, losses in strategies.items():
    v, e = var_es(losses, alpha=0.95)
    print(f"{name:<30} ${losses.mean():>10,.0f} ${v:>10,.0f} ${e:>10,.0f}")

All three strategies spend the SAME $100K; figures below are residual breach loss.
Strategy                            E[Loss]        VaR95         ES95
────────────────────────────── ──────────── ──────────── ────────────
Baseline (no investment)       $   238,630 $ 1,200,000 $ 1,435,928
(a) Spread audit ($25K x 4)    $   192,730 $ 1,200,000 $ 1,358,483
(b) Deep audit ($100K on D)    $   156,430 $ 1,200,000 $ 1,301,996
(c) Contracts ($100K)          $   119,315 $   600,000 $   717,964

# --- Bar chart: residual expected loss and VaR95 by strategy ---

fig, axes = plt.subplots(1, 2, figsize=(10, 5))

strat_labels = ["Baseline", "(a) Spread\naudit", "(b) Deep\naudit", "(c) Contract\ncontrols"]
strat_colors = [LIGHT_GRAY, DARK_BG, DARK_BG, DARK_BG]

ax = axes[0]
means = [l.mean() for l in strategies.values()]
bars = ax.bar(strat_labels, means, color=strat_colors, edgecolor="white")
for bar, m in zip(bars, means):
    ax.text(bar.get_x() + bar.get_width() / 2, m + 1000,
            f"${m:,.0f}", ha="center", va="bottom", fontsize=8, fontweight="bold")
ax.set_ylabel("Residual expected loss")
ax.set_title("Residual Expected Loss")

ax = axes[1]
var95s = [var_es(l, alpha=0.95)[0] for l in strategies.values()]
bars = ax.bar(strat_labels, var95s, color=strat_colors, edgecolor="white")
best_idx = np.argmin(var95s[1:]) + 1  # best VaR95 among the strategies, excluding baseline
bars[best_idx].set_color(ACCENT)
for bar, v in zip(bars, var95s):
    ax.text(bar.get_x() + bar.get_width() / 2, v + 5000,
            f"${v:,.0f}", ha="center", va="bottom", fontsize=8, fontweight="bold")
ax.set_ylabel("Residual VaR at 95%")
ax.set_title("Residual Value-at-Risk (95th percentile)")

for ax in axes:
    ax.yaxis.set_major_formatter(usd_k)

fig.suptitle("Same $100K outlay, three strategies: residual loss after the spend",
             fontsize=13, fontweight="bold", y=1.02)
fig.tight_layout()
plt.show()

print("Bars show breach loss LEFT after the same $100K spend; lower is better.")
print("Under the assumed 50% recovery, contractual controls (c) cut residual expected loss and")
print("VaR95 the most. Deep audit (b) trims ES95 but leaves VaR95 unchanged here. Spread audit")
print("(a) shaves a little off everywhere and barely moves the tail. The ranking is specific to")
print("these assumptions, loss sizes and probabilities -- re-run it for your own vendors.")

Bars show breach loss LEFT after the same $100K spend; lower is better.
Under the assumed 50% recovery, contractual controls (c) cut residual expected loss and
VaR95 the most. Deep audit (b) trims ES95 but leaves VaR95 unchanged here. Spread audit
(a) shaves a little off everywhere and barely moves the tail. The ranking is specific to
these assumptions, loss sizes and probabilities -- re-run it for your own vendors.

6. Real Data: One Application’s Dependency Graph

Sections 1 and 3 used invented vendors and random networks. A software supply chain is a network we can download. deps.dev publishes the resolved dependency graph of open-source packages: every package an install pulls in, and which package requires which.

We take npm react-scripts 5.0.1 – the build toolchain behind many React applications – as one application’s software supply chain. The snapshot in data/ records the graph as deps.dev resolved it on the snapshot date; a fresh install can resolve to different versions. load_dependency_graph(..., refresh=True) downloads the current graph.

from _data import load_dependency_graph

nodes, edges = load_dependency_graph("npm", "react-scripts", "5.0.1")
deps_all = nodes[nodes["relation"] != "SELF"]
# The graph lists some name@version more than once; count each version once.
n_direct = len(deps_all[deps_all["relation"] == "DIRECT"].drop_duplicates(["name", "version"]))
n_all = len(deps_all.drop_duplicates(["name", "version"]))

print(f"Packages the maintainers chose (direct):       {n_direct:,}")
print(f"Distinct package versions an install pulls in: {n_all:,}   ({n_all / n_direct:.0f}x)")
print(f"Graph nodes (some versions listed repeatedly): {len(deps_all):,}")
print(f"Dependency edges:                              {len(edges):,}")

# Section 1's compound exposure, with an illustrative per-package probability
# of a malicious or compromised release in a year.
print("\nP(at least one compromised package), independent packages:")
print(f"{'p per package':>14} {'direct only':>12} {'full graph':>11}")
for p in [1e-5, 1e-4, 1e-3]:
    print(f"{p:>14.0e} {1 - (1 - p) ** n_direct:>12.2%} {1 - (1 - p) ** n_all:>11.1%}")

Packages the maintainers chose (direct):       48
Distinct package versions an install pulls in: 1,210   (25x)
Graph nodes (some versions listed repeatedly): 1,474
Dependency edges:                              2,936

P(at least one compromised package), independent packages:
 p per package  direct only  full graph
         1e-05        0.05%        1.2%
         1e-04        0.48%       11.4%
         1e-03        4.69%       70.2%

Assessing only the packages the maintainers chose understates exposure. For small per-package probabilities the ratio approaches the fan-out; at higher probabilities it shrinks as the full-graph column nears 100%. The per-package probabilities are illustrative.

Blast radius. If one package is compromised (a malicious release), every package that depends on it, directly or through others, pulls it in. A package’s blast radius here is the number of packages in this graph that transitively depend on it – how many other packages one compromise sits underneath.

# Merge repeated copies of a version into one package, and leave out the root.
pkg = nodes["name"] + "@" + nodes["version"]
root = pkg[nodes["relation"] == "SELF"].iloc[0]
links = pd.DataFrame({"child": pkg[edges["to_node"]].values,
                      "parent": pkg[edges["from_node"]].values})
links = links[links["parent"] != root]
parents = links.groupby("child")["parent"].apply(set).to_dict()


def blast_radius(p):
    reached, stack = set(), [p]
    while stack:
        for parent in parents.get(stack.pop(), ()):
            if parent not in reached:
                reached.add(parent)
                stack.append(parent)
    reached.discard(p)  # a package in a dependency cycle does not count itself
    return len(reached)


deps = (nodes[nodes["relation"] != "SELF"].assign(pkg=pkg)
        .sort_values("relation")  # DIRECT sorts first, so a version chosen directly keeps that label
        .drop_duplicates("pkg"))
deps["blast_radius"] = [blast_radius(p) for p in deps["pkg"]]
br = deps["blast_radius"].values

print(f"Blast radius over {len(deps):,} distinct package versions:")
print(f"  Median: {np.median(br):.0f}   P95: {np.percentile(br, 95):.0f}   Max: {br.max()} "
      f"({br.max() / len(deps):.0%} of the graph)")
print(f"  Packages reaching >= 50 others: {(br >= 50).sum()}")
print("\nWidest reach:")
print(deps.sort_values("blast_radius", ascending=False)
      .head(10)[["name", "version", "relation", "blast_radius"]].to_string(index=False))

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.hist(br, bins=np.arange(0, br.max() + 5, 5), color=DARK_BG, edgecolor="white")
ax.axvline(np.median(br), color=LIGHT_GRAY, ls=":", lw=1.5, label=f"Median = {np.median(br):.0f}")
ax.axvline(np.percentile(br, 95), color=ACCENT, ls="--", lw=1.5,
           label=f"95th percentile = {np.percentile(br, 95):.0f}")
ax.set_yscale("log")
ax.set_xlabel("Blast radius (packages that depend on the compromised one)")
ax.set_ylabel("Number of packages (log scale)")
ax.set_title("react-scripts 5.0.1: blast radius of a single compromised package")
ax.legend()
plt.tight_layout()
plt.show()

Blast radius over 1,210 distinct package versions:
  Median: 5   P95: 49   Max: 133 (11% of the graph)
  Packages reaching >= 50 others: 59

Widest reach:
                              name version relation  blast_radius
                        picocolors   1.1.1 INDIRECT           133
                         es-errors   1.3.0 INDIRECT           129
        @babel/helper-plugin-utils  7.29.7 INDIRECT           124
                     function-bind   1.1.2 INDIRECT           124
@babel/helper-validator-identifier  7.29.7 INDIRECT           120
                            hasown   2.0.4 INDIRECT           118
                              gopd   1.2.0 INDIRECT           105
                         js-tokens   4.0.0 INDIRECT           102
                                ms   2.1.3 INDIRECT           102
                es-define-property   1.0.1 INDIRECT           102

Unlike the small random networks in Section 3, where nodes look alike, the real graph is strongly skewed. Most packages sit near the leaves and reach a handful of others; a few small utility packages – none of them chosen directly by the maintainers – sit under up to about a tenth of the tree. As in Section 5, effort pays off most where exposure concentrates: the packages with the widest reach are where review, pinning and provenance checks matter most, and they are packages nobody picked.

Caveats. This is one package’s install, not a whole organization’s software estate, and it counts reach, not damage: a compromised utility deep in a build toolchain may never run in production. The graph is resolved as of the snapshot date.

Data: deps.dev (CC BY 4.0), https://deps.dev.

7. Pitfalls

Independence assumptions dramatically understate supply chain tail risk. Vendors share cloud providers, software dependencies, and threat actor attention. Treating vendor failures as independent is not conservative – it is optimistic. The common-cause model shows how a 15% systemic factor inflates VaR roughly two-fold and creates a fat tail the independence model cannot see.

Vendor count is a risk factor, not just a management challenge. Each additional vendor compounds exposure – sublinearly but relentlessly. Twenty vendors at 5% each produce a 64% compound exposure – not 100%, but far more than the 5% each one contributes individually, and the number only ever climbs. In this model, where any vendor’s breach reaches A, fewer vendors means fewer ways in; consolidation trades that against concentration on the vendors that remain.

Systemic risk cannot be diversified away – it is the opposite of diversification. In financial portfolios, diversification reduces variance. In supply chains, adding vendors adds exposed links: more vendors means more shared infrastructure, more common dependencies, and higher probability that a single systemic event (Log4j, SolarWinds, MOVEit) hits multiple links simultaneously.

Contractual controls (SLAs) shift losses, they do not prevent them. An SLA that guarantees breach notification within 72 hours does not prevent the breach. An indemnification clause can transfer part of the financial loss, up to the vendor’s liability cap and if the vendor can pay; it does not protect your customers’ data. Contracts are loss transfer mechanisms, and treating them as loss prevention overstates their value.

The weakest vendor sets your floor – until you remediate it. If your riskiest vendor has a 12% annual breach probability, your supply chain exposure is anchored near 12% as long as that vendor stays untouched – no amount of hardening the other vendors moves it. But the floor is not fixed: under the assumed 60% reduction, strategy b would take vendor D from 12% to 4.8%; whether targeting the weakest vendor beats spreading the budget depends on audit effects nobody has measured (Section 4) – and on correctly identifying which vendor is weakest.

Cascade models show that the average case hides the tail – plan for the tail. The cascade-size distribution is right-skewed: most cascades fizzle near the seed, but a substantial minority escalate. The expected case is a near-miss; the bad case is most of the network. Risk management is about the bad case.

Most of the software supply chain is indirect. In the react-scripts graph, a few dozen chosen packages pull in more than a thousand others, and the packages with the widest reach are small utilities nobody selected. Assessing only direct dependencies – or only named vendors – misses where the exposure sits.