Operational risk

Interactive exploration of loss distributions and the power law in operational risk

Operational risk is the risk of loss resulting from inadequate or failed internal processes, people, and systems, or from external events (see Hull 2023, chap. 20). Unlike market or credit risk, which a bank consciously assumes, operational risk is a necessary by-product of doing business. Quantifying it requires modeling both how often losses occur and how large they are.

Loss frequency and loss severity

A standard approach to estimating potential operational risk losses combines two distributions:

Loss frequency — How many loss events occur per year? A Poisson distribution with parameter $\lambda$ is commonly used, where $\lambda$ is the expected number of events per year.
Loss severity — Given that a loss occurs, how large is it? A lognormal distribution is often used, parameterized by the mean $\mu$ and standard deviation $\sigma$ of the log-loss.

These two distributions are combined via Monte Carlo simulation. On each trial: (1) sample the number of events $n$ from the frequency distribution; (2) sample $n$ individual losses from the severity distribution; (3) sum them to get the total annual loss. After many trials, we obtain the total loss distribution and can compute risk measures such as the 99.9th percentile (used by regulators for capital).

Note

The Poisson distribution

The Poisson distribution is a natural model for counting rare, independent events over a fixed period. It assumes that in any short time interval there is a small, constant probability of a loss event occurring, and that events are independent of one another. It has a single parameter $\lambda$, which equals both the mean and the variance of the number of events. For example, if $\lambda = 3$, we expect 3 events per year on average, but in any given year we might observe 0, 1, 2, … events with probabilities given by $P(n) = e^{-\lambda}\lambda^n / n!$. If the observed variance of loss counts exceeds the mean, a negative binomial distribution may be more appropriate; if it is less than the mean, a binomial distribution may be used instead.

Interactive simulation

Adjust the parameters below to explore how the loss frequency and loss severity shape the total loss distribution. Observe how heavier tails or higher event frequencies shift and widen the distribution.

Tip

How to experiment

Change $\lambda$ to increase or decrease the average number of events per year.
Adjust $\mu$ and $\sigma$ to shift or widen the severity of individual losses.
Watch how the total loss distribution, its mean, standard deviation, and 99.9th percentile respond.

viewof lambda = Inputs.range([0.5, 10], {
  label: "Loss frequency λ (events/year)",
  step: 0.5,
  value: 3
})

viewof mu = Inputs.range([-1, 2], {
  label: "Log-loss mean μ",
  step: 0.1,
  value: 0
})

viewof sigma = Inputs.range([0.1, 1.5], {
  label: "Log-loss std. dev. σ",
  step: 0.05,
  value: 0.4
})

viewof nYears = Inputs.range([1, 10], {
  label: "Horizon (years)",
  step: 1,
  value: 1
})

viewof nSim = Inputs.range([1000, 100000], {
  label: "Monte Carlo iterations",
  step: 1000,
  value: 20000
})

mutable resimClicks = 0

viewof resimButton = {
  const btn = html`<button style="background: #2f71d5; color: white; border: none; border-radius: 6px; padding: 0.35rem 1rem; font-weight: 500; font-size: 0.9rem; cursor: pointer;">New simulation</button>`
  btn.onclick = () => { mutable resimClicks++ }
  return btn
}

mulberry32 = seed => {
  return function () {
    let t = seed += 0x6D2B79F5
    t = Math.imul(t ^ t >>> 15, t | 1)
    t ^= t + Math.imul(t ^ t >>> 7, t | 61)
    return ((t ^ t >>> 14) >>> 0) / 4294967296
  }
}

randomNormal = rng => {
  let u = 0, v = 0
  while (u === 0) u = rng()
  while (v === 0) v = rng()
  return Math.sqrt(-2.0 * Math.log(u)) * Math.cos(2.0 * Math.PI * v)
}

simSeed = 20260212 + resimClicks

mcResults = {
  const rng = mulberry32(simSeed)
  const totalLosses = []

  for (let i = 0; i < nSim; i++) {
    // Each trial covers nYears years
    let total = 0
    for (let y = 0; y < nYears; y++) {
      // Sample from Poisson(lambda) using inverse transform
      let nEvents = 0
      let p = Math.exp(-lambda)
      let cumP = p
      const u = rng()
      while (u > cumP) {
        nEvents++
        p *= lambda / nEvents
        cumP += p
      }

      // Sample nEvents losses from lognormal(mu, sigma)
      for (let j = 0; j < nEvents; j++) {
        const z = randomNormal(rng)
        total += Math.exp(mu + sigma * z)
      }
    }
    totalLosses.push(total)
  }

  totalLosses.sort((a, b) => a - b)
  return totalLosses
}

mcMean = mcResults.reduce((a, b) => a + b, 0) / mcResults.length

mcSD = {
  const m = mcMean
  const v = mcResults.reduce((a, b) => a + (b - m) ** 2, 0) / (mcResults.length - 1)
  return Math.sqrt(v)
}

mcP999 = mcResults[Math.floor(0.999 * mcResults.length)]
mcP99 = mcResults[Math.floor(0.99 * mcResults.length)]
mcP95 = mcResults[Math.floor(0.95 * mcResults.length)]

// Analytical values for compound Poisson-lognormal over nYears
analyticMean = nYears * lambda * Math.exp(mu + sigma ** 2 / 2)
analyticSD = Math.sqrt(nYears * lambda * Math.exp(2 * mu + 2 * sigma ** 2))

fmtNum = (v, d = 2) => v.toFixed(d)

// --- Total loss distribution histogram ---
lossHist = {
  // Extend slightly beyond the 99.9th percentile so the marker falls within visible bars
  const xMax = mcResults[Math.floor(0.9995 * mcResults.length)]

  const nBins = 60
  const binWidth = xMax / nBins
  const bins = Array.from({ length: nBins }, (_, i) => ({
    x0: i * binWidth,
    x1: (i + 1) * binWidth,
    count: 0
  }))

  for (const v of mcResults) {
    const idx = Math.min(Math.floor(v / binWidth), nBins - 1)
    if (idx >= 0) bins[idx].count++
  }

  // Convert to density
  const totalCount = mcResults.length
  for (const b of bins) {
    b.density = b.count / (totalCount * binWidth)
  }

  return bins
}

severityPDF = {
  const E_L = Math.exp(mu + sigma ** 2 / 2)
  const xMax = Math.min(E_L * 5, Math.exp(mu + 3.5 * sigma))
  const n = 200
  const dx = xMax / n
  const pts = []
  for (let i = 1; i <= n; i++) {
    const x = i * dx
    const logx = Math.log(x)
    const pdf = (1 / (x * sigma * Math.sqrt(2 * Math.PI))) *
      Math.exp(-0.5 * ((logx - mu) / sigma) ** 2)
    pts.push({ x, pdf })
  }
  return pts
}

E_L = Math.exp(mu + sigma ** 2 / 2)
Var_L = (Math.exp(sigma ** 2) - 1) * Math.exp(2 * mu + sigma ** 2)

Plot.plot({
  height: 320,
  marginLeft: 50,
  marginRight: 20,
  x: {
    label: `Total loss over ${nYears} year${nYears > 1 ? "s" : ""} ($ millions)`,
    grid: false
  },
  y: {
    label: "Density",
    grid: true
  },
  marks: [
    Plot.rectY(lossHist, {
      x1: "x0",
      x2: "x1",
      y: "density",
      fill: "#4682b4",
      fillOpacity: 0.7,
      tip: true,
      title: d => `Loss: $${d.x0.toFixed(1)}–${d.x1.toFixed(1)}M\nDensity: ${d.density.toFixed(4)}`
    }),
    Plot.ruleX([mcMean], {
      stroke: "#d62728",
      strokeWidth: 2,
      strokeDash: [6, 3]
    }),
    Plot.ruleX([mcP999], {
      stroke: "#ff7f0e",
      strokeWidth: 2,
      strokeDash: [4, 4]
    }),
    Plot.text([{ x: mcMean, label: "Mean" }], {
      x: "x",
      y: lossHist.reduce((m, b) => Math.max(m, b.density), 0) * 0.95,
      text: "label",
      fill: "#d62728",
      fontWeight: "bold",
      fontSize: 11,
      dx: -20
    }),
    Plot.text([{ x: mcP999, label: "99.9%" }], {
      x: "x",
      y: lossHist.reduce((m, b) => Math.max(m, b.density), 0) * 0.85,
      text: "label",
      fill: "#ff7f0e",
      fontWeight: "bold",
      fontSize: 11,
      dx: -22
    })
  ]
})

html`<table class="table" style="width: 100%;">
<thead><tr><th>Statistic (over ${nYears} year${nYears > 1 ? "s" : ""}, ${nSim.toLocaleString()} iterations)</th><th>Simulation</th><th>Analytical</th></tr></thead>
<tbody>
<tr><td style="font-weight:500;">Mean total loss ($ M)</td><td>${fmtNum(mcMean)}</td><td>${fmtNum(analyticMean)}</td></tr>
<tr><td style="font-weight:500;">Std. dev. total loss ($ M)</td><td>${fmtNum(mcSD)}</td><td>${fmtNum(analyticSD)}</td></tr>
<tr><td style="font-weight:500;">95th percentile ($ M)</td><td>${fmtNum(mcP95)}</td><td style="color:#888;">no closed form</td></tr>
<tr><td style="font-weight:500;">99th percentile ($ M)</td><td>${fmtNum(mcP99)}</td><td style="color:#888;">no closed form</td></tr>
<tr><td style="font-weight:500;">99.9th percentile ($ M)</td><td>${fmtNum(mcP999)}</td><td style="color:#888;">no closed form</td></tr>
</tbody></table>`

The compound Poisson-lognormal distribution has closed-form expressions for its mean and variance, but not for its percentiles. This is precisely why Monte Carlo simulation is needed: it is the most straightforward way to estimate tail risk measures such as the 99.9th percentile that regulators require for capital calculations.

The individual loss severity (not the aggregate) for the current $\mu$ and $\sigma$. The lognormal density has a characteristic right skew — most losses are small, but occasionally a very large loss can occur.

Plot.plot({
  height: 260,
  marginLeft: 50,
  marginRight: 20,
  x: {
    label: "Individual loss ($ millions)",
    grid: false
  },
  y: {
    label: "Density",
    grid: true
  },
  marks: [
    Plot.areaY(severityPDF, {
      x: "x",
      y: "pdf",
      fill: "#4682b4",
      fillOpacity: 0.2,
      curve: "natural"
    }),
    Plot.line(severityPDF, {
      x: "x",
      y: "pdf",
      stroke: "#4682b4",
      strokeWidth: 2,
      curve: "natural"
    })
  ]
})

html`<table class="table" style="width: 100%;">
<thead><tr><th colspan="2">Individual loss severity (lognormal)</th></tr></thead>
<tbody>
<tr><td style="font-weight:500;">Expected individual loss ($ M)</td><td>${fmtNum(E_L, 3)}</td></tr>
<tr><td style="font-weight:500;">Std. dev. of individual loss ($ M)</td><td>${fmtNum(Math.sqrt(Var_L), 3)}</td></tr>
</tbody></table>`

The power law and heavy tails

For large losses, operational risk data often follow a power law:

\[ \text{Prob}(v > x) = K x^{-\alpha} \]

where $v$ is the loss, $x$ is a threshold, and $K$ and $\alpha$ are constants. The exponent $\alpha$ determines the heaviness of the tail: a smaller $\alpha$ means a heavier tail and dramatically larger extreme losses.

A key practical implication is that if we know the probability of exceeding one threshold, we can estimate the probability of exceeding any other threshold using the ratio property:

\[ \text{Prob}(v > x_2) = \text{Prob}(v > x_1) \times \left(\frac{x_2}{x_1}\right)^{-\alpha} \]

Interactive power law explorer

Set a known exceedance probability and threshold below, then vary $\alpha$ to see how the tail probabilities and extreme quantiles change. Notice how sensitive the results are to $\alpha$.

Tip

How to experiment

Set the baseline: “There is a probability $p$ of losses exceeding $x_1$.”
Drag $\alpha$ to see how the tail becomes heavier or lighter.
Read off the probability of exceeding multiples of the baseline, and the loss at key confidence levels (99%, 99.9%, 99.97%).

viewof x1 = Inputs.range([1, 100], {
  label: "Baseline threshold x₁ ($ millions)",
  step: 1,
  value: 20
})

viewof p1 = Inputs.range([0.01, 0.30], {
  label: "Prob(loss > x₁)",
  step: 0.01,
  value: 0.10
})

viewof alpha = Inputs.range([0.2, 2.0], {
  label: "Power law exponent α",
  step: 0.05,
  value: 0.8
})

// Exceedance probabilities at multiples of x1
plMultiples = {
  const multiples = [1, 2, 4, 5, 10, 20, 50]
  return multiples.map(m => {
    const x = x1 * m
    const prob = p1 * Math.pow(m, -alpha)
    return { multiple: m, x, prob: Math.min(prob, 1) }
  })
}

quantile_from_power = target => {
  // P(v > x) = p1 * (x/x1)^{-alpha} = target
  // (x/x1)^alpha = p1/target
  // x = x1 * (p1/target)^{1/alpha}
  return x1 * Math.pow(p1 / target, 1 / alpha)
}

q99 = quantile_from_power(0.01)
q999 = quantile_from_power(0.001)
q9997 = quantile_from_power(0.0003)

powerLawCurve = {
  const pts = []
  const xMin = x1
  const xMax = quantile_from_power(0.0001)
  const cappedMax = Math.min(xMax, x1 * 200)
  const n = 300
  for (let i = 0; i <= n; i++) {
    // Use log spacing for better coverage of the tail
    const logX = Math.log(xMin) + (Math.log(cappedMax) - Math.log(xMin)) * (i / n)
    const x = Math.exp(logX)
    const prob = p1 * Math.pow(x / x1, -alpha)
    pts.push({ x, prob: Math.min(prob, 1) })
  }
  return pts
}

alphaComparison = {
  const alphas = [0.3, 0.5, 0.8, 1.0, 1.5]
  const pts = []
  const xMin = x1
  const xMax = x1 * 100
  const n = 200
  for (const a of alphas) {
    for (let i = 0; i <= n; i++) {
      const logX = Math.log(xMin) + (Math.log(xMax) - Math.log(xMin)) * (i / n)
      const x = Math.exp(logX)
      const prob = p1 * Math.pow(x / x1, -a)
      if (prob >= 1e-6 && prob <= 1) {
        pts.push({ x, prob, alpha: `α = ${a}` })
      }
    }
  }
  return pts
}

fmtPct = v => {
  const pct = v * 100
  if (pct >= 1) return `${pct.toFixed(0)}%`
  if (pct >= 0.01) return `${pct.toPrecision(2)}%`
  return `${pct.toExponential(0)}%`
}

Plot.plot({
  height: 320,
  marginLeft: 60,
  marginRight: 20,
  x: {
    label: "Loss threshold x ($ millions)",
    type: "log"
  },
  y: {
    label: "P(loss > x)",
    type: "log",
    grid: true,
    tickFormat: fmtPct
  },
  marks: [
    Plot.line(powerLawCurve, {
      x: "x",
      y: "prob",
      stroke: "#d62728",
      strokeWidth: 2.5
    }),
    Plot.dot(plMultiples.filter(d => d.prob > 0.00005 && d.prob <= 1), {
      x: "x",
      y: "prob",
      fill: "#d62728",
      r: 4,
      tip: true,
      title: d => `${d.multiple}× threshold: $${d.x}M\nP(loss > $${d.x}M) = ${(d.prob * 100).toFixed(4)}%`
    }),
    Plot.ruleY([0.01], { stroke: "#888", strokeDash: [3, 3], strokeOpacity: 0.5 }),
    Plot.ruleY([0.001], { stroke: "#888", strokeDash: [3, 3], strokeOpacity: 0.5 }),
    Plot.ruleY([0.0003], { stroke: "#888", strokeDash: [3, 3], strokeOpacity: 0.5 }),
    Plot.text([{ y: 0.01, label: "1%" }], {
      y: "y",
      x: powerLawCurve[powerLawCurve.length - 1].x,
      text: "label",
      fill: "#888",
      fontSize: 10,
      dx: 5,
      textAnchor: "start"
    }),
    Plot.text([{ y: 0.001, label: "0.1%" }], {
      y: "y",
      x: powerLawCurve[powerLawCurve.length - 1].x,
      text: "label",
      fill: "#888",
      fontSize: 10,
      dx: 5,
      textAnchor: "start"
    }),
    Plot.text([{ y: 0.0003, label: "0.03%" }], {
      y: "y",
      x: powerLawCurve[powerLawCurve.length - 1].x,
      text: "label",
      fill: "#888",
      fontSize: 10,
      dx: 5,
      textAnchor: "start"
    })
  ]
})

html`<table class="table" style="width: 100%;">
<thead><tr>
  <th>Multiple</th>
  <th>Loss threshold ($ M)</th>
  <th>P(loss > threshold)</th>
</tr></thead>
<tbody>
${plMultiples.filter(d => d.prob <= 1 && d.prob >= 0.00001).map(d =>
  `<tr>
    <td>${d.multiple}×</td>
    <td style="text-align:right;">$${d.x.toFixed(0)}</td>
    <td style="text-align:right;">${(d.prob * 100).toFixed(4)}%</td>
  </tr>`
).join("")}
</tbody></table>`

The loss levels at key confidence levels — the amounts exceeded with a given (small) probability.

html`<table class="table" style="width: 100%;">
<thead><tr>
  <th>Confidence level</th>
  <th>Exceedance probability</th>
  <th>Loss ($ millions)</th>
</tr></thead>
<tbody>
<tr><td>99%</td><td>1%</td><td style="text-align:right;">${fmtNum(q99, 1)}</td></tr>
<tr><td>99.9%</td><td>0.1%</td><td style="text-align:right;">${fmtNum(q999, 1)}</td></tr>
<tr><td>99.97%</td><td>0.03%</td><td style="text-align:right;">${fmtNum(q9997, 1)}</td></tr>
</tbody></table>`

Exceedance probability curves for several values of $\alpha$, holding the baseline fixed. A lower $\alpha$ produces a dramatically heavier tail.

Plot.plot({
  height: 340,
  marginLeft: 60,
  marginRight: 20,
  color: {
    legend: true
  },
  x: {
    label: "Loss threshold x ($ millions)",
    type: "log"
  },
  y: {
    label: "P(loss > x)",
    type: "log",
    grid: true,
    tickFormat: fmtPct
  },
  marks: [
    Plot.line(alphaComparison, {
      x: "x",
      y: "prob",
      stroke: "alpha",
      strokeWidth: 2
    }),
    Plot.ruleY([0.001], { stroke: "#888", strokeDash: [3, 3], strokeOpacity: 0.4 }),
    Plot.text([{ y: 0.001, label: "0.1% (99.9th percentile)" }], {
      y: "y",
      x: alphaComparison[alphaComparison.length - 1].x * 0.8,
      text: "label",
      fill: "#888",
      fontSize: 10,
      textAnchor: "end"
    })
  ]
})

Note

Why $\alpha$ matters so much

With $\alpha = 0.5$ and a 10% chance of exceeding $20M, the 99.9th percentile loss is in the billions. With $\alpha = 1.5$, it drops to under $100M. Accurately estimating $\alpha$ from data is critical for operational risk capital calculations. The power law implies that the probability of extreme losses declines polynomially — much more slowly than the exponential decay of normal or lognormal tails.

References

Hull, John. 2023. Risk Management and Financial Institutions. 6th ed. New Jersey: John Wiley & Sons.