Filtered Historical Simulation

Interactive exploration of Filtered Historical Simulation (FHS), combining GARCH volatility dynamics with the empirical distribution of standardized shocks

Filtered Historical Simulation (FHS) is a semi-parametric method for computing Value-at-Risk and Expected Shortfall. It combines the best of two worlds: GARCH-type models for capturing time-varying volatility, and the empirical distribution of standardized shocks for capturing fat tails and skewness without distributional assumptions (Christoffersen 2012, chap. 6; Barone-Adesi et al. 1999).

Recall that GARCH models decompose daily returns as

\[ R_{PF,t+1} = \sigma_{PF,t+1}\, z_{t+1}, \quad z_{t+1} \sim \text{i.i.d.}\;D(0,1) \]

where \(\sigma_{PF,t+1}\) is the conditional standard deviation and \(z_{t+1}\) is the standardized innovation. Parametric approaches assume a specific distribution \(D\) (e.g., normal or Student-t), while standard Historical Simulation ignores the time-varying volatility altogether. FHS takes a different path:

Estimate a GARCH model to obtain the conditional volatility series \(\{\hat{\sigma}_t\}\)
Filter the returns to extract standardized residuals \(\hat{z}_t = R_t / \hat{\sigma}_t\), which should be approximately i.i.d.
Compute risk measures using the empirical distribution of the \(\hat{z}_t\), scaled by the current volatility forecast \(\hat{\sigma}_{t+1}\)

The FHS VaR at tail probability \(p\) is

\[ \text{VaR}_{t+1}^{p} = -\hat{\sigma}_{PF,t+1} \cdot \text{Percentile}\left(\{\hat{z}_{t+1-\tau}\}_{\tau=1}^{m},\; 100p\right) \]

and the Expected Shortfall is

\[ \text{ES}_{t+1}^{p} = -\hat{\sigma}_{PF,t+1} \cdot \frac{1}{|\{i : \hat{z}_i < \hat{z}_{(p)}\}|} \sum_{\hat{z}_i < \hat{z}_{(p)}} \hat{z}_i \]

Note

A semi-parametric approach. FHS is parametric for the volatility dynamics (GARCH) but nonparametric for the shock distribution (empirical). This means it inherits the volatility responsiveness of GARCH and the distributional flexibility of Historical Simulation.

// ============================================================
// MATH UTILITIES: Normal distribution
// ============================================================

// Standard normal CDF (Abramowitz & Stegun approximation)
normalCDF = x => {
  const a1 = 0.254829592, a2 = -0.284496736, a3 = 1.421413741
  const a4 = -1.453152027, a5 = 1.061405429, p = 0.3275911
  const sign = x < 0 ? -1 : 1
  const z = Math.abs(x) / Math.sqrt(2)
  const t = 1.0 / (1.0 + p * z)
  const y = 1 - (((((a5 * t + a4) * t) + a3) * t + a2) * t + a1) * t * Math.exp(-z * z)
  return 0.5 * (1 + sign * y)
}

// Standard normal PDF
normalPDF = x => Math.exp(-x * x / 2) / Math.sqrt(2 * Math.PI)

// Inverse normal CDF (Acklam's algorithm, max |error| ~ 1.15e-9)
qnorm = {
  const a1 = -3.969683028665376e+01, a2 =  2.209460984245205e+02
  const a3 = -2.759285104469687e+02, a4 =  1.383577518672690e+02
  const a5 = -3.066479806614716e+01, a6 =  2.506628277459239e+00
  const b1 = -5.447609879822406e+01, b2 =  1.615858368580409e+02
  const b3 = -1.556989798598866e+02, b4 =  6.680131188771972e+01
  const b5 = -1.328068155288572e+01
  const c1 = -7.784894002430293e-03, c2 = -3.223964580411365e-01
  const c3 = -2.400758277161838e+00, c4 = -2.549732539343734e+00
  const c5 =  4.374664141464968e+00, c6 =  2.938163982698783e+00
  const d1 =  7.784695709041462e-03, d2 =  3.224671290700398e-01
  const d3 =  2.445134137142996e+00, d4 =  3.754408661907416e+00
  const pLow = 0.02425, pHigh = 1 - pLow
  return p => {
    if (p <= 0) return -Infinity
    if (p >= 1) return Infinity
    if (p < pLow) {
      const q = Math.sqrt(-2 * Math.log(p))
      return (((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1)
    }
    if (p <= pHigh) {
      const q = p - 0.5, r = q * q
      return (((((a1*r+a2)*r+a3)*r+a4)*r+a5)*r+a6)*q / (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r+1)
    }
    const q = Math.sqrt(-2 * Math.log(1 - p))
    return -(((((c1*q+c2)*q+c3)*q+c4)*q+c5)*q+c6) / ((((d1*q+d2)*q+d3)*q+d4)*q+1)
  }
}

// Normal ES (as a positive number): ES = σ × φ(Φ⁻¹(p)) / p
normalES = (p, sigma) => sigma * normalPDF(qnorm(p)) / p

// ============================================================
// MATH UTILITIES: Log-gamma and Student-t distribution
// ============================================================

// Log-gamma function (Lanczos approximation, g=7)
lgamma = {
  const g = 7
  const coef = [
    0.99999999999980993, 676.5203681218851, -1259.1392167224028,
    771.32342877765313, -176.61502916214059, 12.507343278686905,
    -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7
  ]
  function _lgamma(x) {
    if (x <= 0) return Infinity
    if (x < 0.5) {
      return Math.log(Math.PI / Math.sin(Math.PI * x)) - _lgamma(1 - x)
    }
    x -= 1
    let a = coef[0]
    const t = x + g + 0.5
    for (let i = 1; i < g + 2; i++) a += coef[i] / (x + i)
    return 0.5 * Math.log(2 * Math.PI) + (x + 0.5) * Math.log(t) - t + Math.log(a)
  }
  return _lgamma
}

// Standard (non-standardized) Student-t PDF
tPDF = (x, d) => {
  return Math.exp(
    lgamma((d + 1) / 2) - lgamma(d / 2)
    - 0.5 * Math.log(d * Math.PI)
    - ((d + 1) / 2) * Math.log(1 + x * x / d)
  )
}

// Regularized incomplete beta function (Numerical Recipes algorithm)
regIncBeta = {
  function betacf(a, b, x) {
    const qab = a + b, qap = a + 1, qam = a - 1
    const TINY = 1e-30, EPS = 1e-14, MAXITER = 200
    let c = 1
    let d = 1 - qab * x / qap
    if (Math.abs(d) < TINY) d = TINY
    d = 1 / d
    let h = d
    for (let m = 1; m <= MAXITER; m++) {
      const m2 = 2 * m
      let aa = m * (b - m) * x / ((qam + m2) * (a + m2))
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY
      d = 1 / d; h *= d * c
      aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2))
      d = 1 + aa * d; if (Math.abs(d) < TINY) d = TINY
      c = 1 + aa / c; if (Math.abs(c) < TINY) c = TINY
      d = 1 / d
      const del = d * c
      h *= del
      if (Math.abs(del - 1) < EPS) break
    }
    return h
  }
  function _rib(a, b, x) {
    if (x < 0 || x > 1) return NaN
    if (x === 0) return 0
    if (x === 1) return 1
    if (x > (a + 1) / (a + b + 2)) {
      return 1 - _rib(b, a, 1 - x)
    }
    const lbeta = lgamma(a) + lgamma(b) - lgamma(a + b)
    const front = Math.exp(Math.log(x) * a + Math.log(1 - x) * b - lbeta) / a
    return front * betacf(a, b, x)
  }
  return _rib
}

// Student-t CDF
tCDF = (x, d) => {
  const xt = d / (d + x * x)
  if (x >= 0) {
    return 1 - 0.5 * regIncBeta(d / 2, 0.5, xt)
  } else {
    return 0.5 * regIncBeta(d / 2, 0.5, xt)
  }
}

// Student-t quantile via Newton-Raphson
qt = (p, d) => {
  if (p <= 0) return -Infinity
  if (p >= 1) return Infinity
  if (p === 0.5) return 0
  let x = qnorm(p)
  for (let i = 0; i < 50; i++) {
    const err = tCDF(x, d) - p
    const deriv = tPDF(x, d)
    if (Math.abs(deriv) < 1e-30) break
    const dx = err / deriv
    x -= dx
    if (Math.abs(dx) < 1e-12 * (1 + Math.abs(x))) break
  }
  return x
}

// Normalizing constant C(d) for standardized t
stdTC = d => Math.exp(lgamma((d + 1) / 2) - lgamma(d / 2) - 0.5 * Math.log(Math.PI * (d - 2)))

// Standardized t PDF: f(z; d) = C(d) * (1 + z²/(d-2))^(-(1+d)/2)
stdTPDF = (z, d) => {
  const C = stdTC(d)
  return C * Math.pow(1 + z * z / (d - 2), -(1 + d) / 2)
}

// Standardized t quantile
stdTQuantile = (p, d) => Math.sqrt((d - 2) / d) * qt(p, d)

// Standardized t CDF
stdTCDF = (z, d) => tCDF(z * Math.sqrt(d / (d - 2)), d)

// ============================================================
// PRNG UTILITIES
// ============================================================

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

  function boxMuller(rng) {
    const u1 = rng(), u2 = rng()
    return Math.sqrt(-2 * Math.log(u1)) * Math.cos(2 * Math.PI * u2)
  }

  function gammaSample(rng, shape) {
    if (shape < 1) {
      return gammaSample(rng, shape + 1) * Math.pow(rng(), 1 / shape)
    }
    const d = shape - 1 / 3
    const c = 1 / Math.sqrt(9 * d)
    while (true) {
      let x, v
      do {
        x = boxMuller(rng)
        v = 1 + c * x
      } while (v <= 0)
      v = v * v * v
      const u = rng()
      if (u < 1 - 0.0331 * (x * x) * (x * x)) return d * v
      if (Math.log(u) < 0.5 * x * x + d * (1 - v + Math.log(v))) return d * v
    }
  }

  function chiSquaredSample(rng, d) {
    return 2 * gammaSample(rng, d / 2)
  }

  function stdTSample(rng, d) {
    const normal = boxMuller(rng)
    const chi2 = chiSquaredSample(rng, d)
    return normal / Math.sqrt(chi2 / d) * Math.sqrt((d - 2) / d)
  }

  return { mulberry32, boxMuller, stdTSample }
}

// ============================================================
// TIME SERIES AND FORMATTING UTILITIES
// ============================================================

// Autocorrelation function (for lags 1..maxLag)
computeACF = (series, maxLag) => {
  const n = series.length
  const mean = series.reduce((a, b) => a + b, 0) / n
  const demeaned = series.map(x => x - mean)
  const var0 = demeaned.reduce((a, b) => a + b * b, 0) / n
  const result = []
  for (let k = 1; k <= maxLag; k++) {
    let sum = 0
    for (let t = k; t < n; t++) sum += demeaned[t] * demeaned[t - k]
    result.push({ lag: k, acf: sum / (n * var0) })
  }
  return result
}

// Ljung-Box statistic
ljungBox = (series, K) => {
  const acfVals = computeACF(series, K)
  const n = series.length
  let stat = 0
  for (const { lag, acf } of acfVals) {
    stat += ((n + 2) / (n - lag)) * acf * acf
  }
  return n * stat
}

// Descriptive statistics helper
descStats = (arr) => {
  const n = arr.length
  const mean = arr.reduce((a, b) => a + b, 0) / n
  const m2 = arr.reduce((a, b) => a + (b - mean) ** 2, 0) / n
  const sd = Math.sqrt(m2)
  const z = arr.map(x => (x - mean) / sd)
  const skew = z.reduce((a, b) => a + b ** 3, 0) / n
  const kurt = z.reduce((a, b) => a + b ** 4, 0) / n
  return { mean, variance: m2, sd, skew, kurt }
}

// Format number to fixed decimals
fmt = (x, d) => x === undefined || isNaN(x) ? "N/A" : x.toFixed(d)

// Percentage formatter for plot axes
pctFmt = x => (x * 100).toFixed(1) + "%"

// Round to 10 decimal places to avoid floating-point artifacts
fpRound = x => Math.round(x * 1e10) / 1e10

// ============================================================
// LEGEND HELPER (interactive: click items to toggle series)
// ============================================================
// Returns a viewof-compatible DOM element. Its .value is a Set
// of hidden keys. Clicking an item toggles it and dispatches an
// "input" event so dependent OJS cells (plots) re-render.
// Each item: { key, label, color, type }. key defaults to label.
legend = (items) => {
  const el = document.createElement("div")
  el.style.cssText = "display:flex; flex-wrap:wrap; margin-top:-4px; margin-bottom:6px;"
  const hidden = new Set()

  for (const d of items) {
    const key = d.key || d.label
    const span = document.createElement("span")
    span.style.cssText = "display:inline-flex; align-items:center; gap:4px; margin-right:14px; cursor:pointer; user-select:none; transition:opacity 0.15s;"

    let swatchHTML
    if (d.type === "dot") {
      swatchHTML = `<svg width="12" height="12"><circle cx="6" cy="6" r="5" fill="${d.color}" opacity="0.8"/></svg>`
    } else if (d.type === "dashed") {
      swatchHTML = `<svg width="22" height="12"><line x1="0" y1="6" x2="22" y2="6" stroke="${d.color}" stroke-width="2" stroke-dasharray="4 2"/></svg>`
    } else if (d.type === "rect") {
      swatchHTML = `<svg width="14" height="14"><rect width="14" height="14" fill="${d.color}"/></svg>`
    } else {
      swatchHTML = `<svg width="22" height="12"><line x1="0" y1="6" x2="22" y2="6" stroke="${d.color}" stroke-width="2"/></svg>`
    }

    span.innerHTML = `${swatchHTML}<span style="font-size:0.82rem;">${d.label}</span>`
    span.addEventListener("click", () => {
      const nowHidden = !hidden.has(key)
      if (nowHidden) hidden.add(key); else hidden.delete(key)
      span.style.opacity = nowHidden ? "0.35" : "1"
      span.querySelector("span").style.textDecoration = nowHidden ? "line-through" : "none"
      el.value = new Set(hidden)
      el.dispatchEvent(new Event("input", {bubbles: true}))
    })
    el.appendChild(span)
  }

  el.value = new Set(hidden)
  return el
}

1. The filtering step

The key idea behind FHS is that while raw returns \(R_t\) exhibit volatility clustering (they are not i.i.d.), the standardized residuals

\[ \hat{z}_t = \frac{R_t}{\hat{\sigma}_t} \]

should be approximately i.i.d. if the GARCH model is correctly specified. This “filtering” step removes the time-varying volatility from the returns, leaving behind the pure shock component whose empirical distribution can be used for simulation.

Note

Using true parameters. In this simulation we use the true GARCH parameters (since we generated the data). In practice, the parameters would be estimated via maximum likelihood (see the MLE page and the GARCH models page).

Tip

How to experiment

Start with Normal shocks and observe that both raw returns and standardized residuals look similar in the QQ plot. Then switch to “Standardized t” and lower the degrees of freedom: the raw returns will show exaggerated fat tails (S-shaped QQ plot) due to the combination of fat-tailed shocks and volatility clustering, while the standardized residuals reveal only the true shock distribution. Check the ACF tab to confirm that filtering removes the autocorrelation from squared returns.

viewof fsAlpha = Inputs.range([0.01, 0.25], {
  label: "α (GARCH reaction)",
  step: 0.01,
  value: 0.10
})

viewof fsBeta = Inputs.range([0.50, 0.98], {
  label: "β (GARCH persistence)",
  step: 0.01,
  value: 0.85
})

viewof fsLRVol = Inputs.range([0.5, 3.0], {
  label: "Long-run daily volatility (%)",
  step: 0.1,
  value: 1.5
})

viewof fsDist = Inputs.radio(["Normal", "Standardized t"], {
  label: "True shock distribution",
  value: "Normal"
})

viewof fsDf = fsDist === "Standardized t" ? Inputs.range([4, 30], {
  label: "Degrees of freedom d",
  step: 1,
  value: 8
}) : Object.assign(html`<span></span>`, {value: 8})

viewof fsN = Inputs.range([500, 3000], {
  label: "Sample size",
  step: 100,
  value: 2000
})

mutable fsSeed = 0

viewof fsReseed = {
  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 sample</button>`
  btn.onclick = () => { mutable fsSeed++ }
  return btn
}

fsPers = fsAlpha + fsBeta
fsVL = (fsLRVol / 100) ** 2
fsOmega = fsVL * (1 - fsAlpha - fsBeta)

// Simulate GARCH(1,1) process
fsSim = {
  fsSeed  // reactivity trigger
  const rng = rng_utils.mulberry32(42 + fsSeed)
  const n = fsN
  const alpha = fsAlpha, beta = fsBeta, omega = fsOmega, vl = fsVL

  const ret = new Array(n)
  const sig2 = new Array(n)
  const shock = new Array(n)

  sig2[0] = vl
  for (let t = 0; t < n; t++) {
    const z = fsDist === "Normal"
      ? rng_utils.boxMuller(rng)
      : rng_utils.stdTSample(rng, fsDf)
    shock[t] = z
    ret[t] = Math.sqrt(sig2[t]) * z
    if (t < n - 1) {
      sig2[t + 1] = omega + alpha * ret[t] * ret[t] + beta * sig2[t]
    }
  }

  // Standardized residuals (using true parameters)
  const stdResid = ret.map((r, i) => r / Math.sqrt(sig2[i]))

  return { ret, sig2, shock, stdResid }
}

// Plot data
fsPlotData = fsSim.ret.map((r, i) => ({
  t: i + 1,
  ret: r,
  vol: Math.sqrt(fsSim.sig2[i]),
  stdResid: fsSim.stdResid[i],
  retSq: r * r,
  stdResidSq: fsSim.stdResid[i] * fsSim.stdResid[i]
}))

fsACFLags = 30
fsACFRetSq = computeACF(fsSim.ret.map(r => r * r), fsACFLags)
fsACFStdSq = computeACF(fsSim.stdResid.map(z => z * z), fsACFLags)
fsBartlett = 2 / Math.sqrt(fsN)

fsLBRetSq = ljungBox(fsSim.ret.map(r => r * r), 15)
fsLBStdSq = ljungBox(fsSim.stdResid.map(z => z * z), 15)

fsStatsRet = descStats(fsSim.ret)
fsStatsStd = descStats(fsSim.stdResid)

{
  const retPlot = Plot.plot({
    height: 200,
    marginLeft: 60,
    marginBottom: 5,
    x: { label: null, grid: false },
    y: { label: "Return", grid: true, tickFormat: pctFmt },
    marks: [
      Plot.areaY(fsPlotData, {
        x: "t",
        y1: d => 2 * d.vol,
        y2: d => -2 * d.vol,
        fill: "#e0e0e0",
        fillOpacity: 0.5
      }),
      Plot.ruleY([0], { stroke: "#ccc" }),
      Plot.line(fsPlotData, { x: "t", y: "ret", stroke: "#333", strokeWidth: 0.5 })
    ]
  })

  const stdPlot = Plot.plot({
    height: 200,
    marginLeft: 60,
    marginTop: 5,
    x: { label: "Day", grid: false },
    y: { label: "Standardized shock", grid: true },
    marks: [
      Plot.ruleY([2, -2], { stroke: "#999", strokeDasharray: "6 3" }),
      Plot.ruleY([0], { stroke: "#ccc" }),
      Plot.line(fsPlotData, { x: "t", y: "stdResid", stroke: "#2f71d5", strokeWidth: 0.5 })
    ]
  })

  return html`<div>
    ${retPlot}
    ${stdPlot}
  </div>`
}

html`<p style="color:#666;font-size:0.85rem;"><strong>Top:</strong> Raw returns R<sub>t</sub> with ±2σ<sub>t</sub> bands (gray shading). Volatility clustering is visible. <strong>Bottom:</strong> Standardized shocks ẑ<sub>t</sub> = R<sub>t</sub>/σ<sub>t</sub> with ±2 reference lines. The shocks look much more homogeneous.</p>`

{
  const retSorted = fsSim.ret.slice().sort((a, b) => a - b)
  const stdSorted = fsSim.stdResid.slice().sort((a, b) => a - b)
  const n = retSorted.length

  // Standardize raw returns for QQ plot
  const retMean = retSorted.reduce((a, b) => a + b, 0) / n
  const retSD = Math.sqrt(retSorted.reduce((a, b) => a + (b - retMean) ** 2, 0) / (n - 1))
  const retStd = retSorted.map(r => (r - retMean) / retSD)

  const qqRet = retStd.map((z, i) => ({
    theoretical: qnorm((i + 0.5) / n),
    empirical: z
  }))

  const qqStd = stdSorted.map((z, i) => ({
    theoretical: qnorm((i + 0.5) / n),
    empirical: z
  }))

  const lim = Math.max(4, Math.max(
    Math.abs(qqRet[0].empirical), Math.abs(qqRet[n - 1].empirical),
    Math.abs(qqStd[0].empirical), Math.abs(qqStd[n - 1].empirical)
  ))

  const retQQ = Plot.plot({
    width: 350, height: 320, marginLeft: 50,
    x: { label: "Normal quantile", domain: [-lim, lim] },
    y: { label: "Empirical quantile", domain: [-lim, lim], grid: true },
    marks: [
      Plot.line([[-lim, -lim], [lim, lim]], { stroke: "#d62728", strokeWidth: 1.5 }),
      Plot.dot(qqRet, { x: "theoretical", y: "empirical", r: 1, fill: "#333", fillOpacity: 0.3 })
    ]
  })

  const stdQQ = Plot.plot({
    width: 350, height: 320, marginLeft: 50,
    x: { label: "Normal quantile", domain: [-lim, lim] },
    y: { label: "Empirical quantile", domain: [-lim, lim], grid: true },
    marks: [
      Plot.line([[-lim, -lim], [lim, lim]], { stroke: "#d62728", strokeWidth: 1.5 }),
      Plot.dot(qqStd, { x: "theoretical", y: "empirical", r: 1, fill: "#2f71d5", fillOpacity: 0.3 })
    ]
  })

  const leftDiv = html`<div style="flex:1;min-width:300px;max-width:400px;"><h4 style="text-align:center;margin-bottom:4px;">Raw returns</h4></div>`
  leftDiv.appendChild(retQQ)
  const rightDiv = html`<div style="flex:1;min-width:300px;max-width:400px;"><h4 style="text-align:center;margin-bottom:4px;">Standardized shocks</h4></div>`
  rightDiv.appendChild(stdQQ)

  const container = html`<div style="display:flex;gap:20px;flex-wrap:wrap;"></div>`
  container.appendChild(leftDiv)
  container.appendChild(rightDiv)
  return container
}

html`<p style="color:#666;font-size:0.85rem;"><strong>Left:</strong> QQ plot of raw returns (standardized to zero mean, unit variance) vs. normal. The S-shape results from mixing different volatility regimes. <strong>Right:</strong> QQ plot of standardized shocks vs. normal. Departures from the 45-degree line now reflect the true shock distribution only.</p>`

viewof fsACFLegend = legend([
  { key: "raw", label: "R² (raw returns)", color: "#333", type: "rect" },
  { key: "std", label: "ẑ² (standardized shocks)", color: "#2f71d5", type: "rect" },
  { key: "bartlett", label: "95% Bartlett band", color: "#d62728", type: "dashed" }
])

{
  const rawBars = fsACFRetSq.map(d => ({ lag: d.lag, acf: d.acf }))
  const stdBars = fsACFStdSq.map(d => ({ lag: d.lag, acf: d.acf }))
  const h = fsACFLegend

  return Plot.plot({
    height: 350, marginLeft: 60,
    x: { label: "Lag", grid: false, ticks: [5, 10, 15, 20, 25, 30] },
    y: { label: "Autocorrelation", grid: true, domain: [-0.1, Math.max(0.3, ...fsACFRetSq.map(d => d.acf)) * 1.1] },
    marks: [
      Plot.ruleY([0], { stroke: "#ccc" }),
      ...(!h.has("bartlett") ? [
        Plot.ruleY([fsBartlett], { stroke: "#d62728", strokeDasharray: "4 3", strokeWidth: 1 }),
        Plot.ruleY([-fsBartlett], { stroke: "#d62728", strokeDasharray: "4 3", strokeWidth: 1 })
      ] : []),
      ...(!h.has("raw") ? [
        Plot.rectY(rawBars, { x1: d => d.lag - 0.4, x2: d => d.lag, y: "acf", fill: "#333" })
      ] : []),
      ...(!h.has("std") ? [
        Plot.rectY(stdBars, { x1: d => d.lag, x2: d => d.lag + 0.4, y: "acf", fill: "#2f71d5" })
      ] : [])
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">Ljung-Box(15) for R²: <strong>${fmt(fsLBRetSq, 1)}</strong> | Ljung-Box(15) for ẑ²: <strong>${fmt(fsLBStdSq, 1)}</strong> | χ²(15) critical value at 5%: 24.996. Strong autocorrelation in squared returns (volatility clustering) is removed after GARCH filtering.</p>`

html`<table class="table" style="width:100%;">
<thead><tr>
  <th>Statistic</th>
  <th>Raw returns R<sub>t</sub></th>
  <th>Standardized shocks ẑ<sub>t</sub></th>
</tr></thead>
<tbody>
<tr><td style="font-weight:500;">Mean</td>
    <td>${fmt(fsStatsRet.mean * 100, 4)}%</td>
    <td>${fmt(fsStatsStd.mean, 4)}</td></tr>
<tr><td style="font-weight:500;">Std. deviation</td>
    <td>${fmt(fsStatsRet.sd * 100, 4)}%</td>
    <td>${fmt(fsStatsStd.sd, 4)}</td></tr>
<tr><td style="font-weight:500;">Skewness</td>
    <td>${fmt(fsStatsRet.skew, 4)}</td>
    <td>${fmt(fsStatsStd.skew, 4)}</td></tr>
<tr><td style="font-weight:500;">Kurtosis</td>
    <td>${fmt(fsStatsRet.kurt, 4)}</td>
    <td>${fmt(fsStatsStd.kurt, 4)}</td></tr>
<tr><td style="font-weight:500;">Excess kurtosis</td>
    <td style="font-weight:700; color:#d62728;">${fmt(fsStatsRet.kurt - 3, 4)}</td>
    <td style="font-weight:700; color:#d62728;">${fmt(fsStatsStd.kurt - 3, 4)}</td></tr>
<tr><td style="font-weight:500;">Ljung-Box(15) of squared series</td>
    <td style="font-weight:700;">${fmt(fsLBRetSq, 1)}</td>
    <td>${fmt(fsLBStdSq, 1)}</td></tr>
</tbody></table>
<p style="color:#666;font-size:0.85rem;">Raw returns have higher excess kurtosis than the standardized shocks because volatility clustering adds additional leptokurtosis beyond what the shock distribution contributes. The Ljung-Box statistic confirms that GARCH filtering removes the autocorrelation from squared returns.</p>`

2. FHS VaR and ES

Once we have extracted the standardized shocks \(\hat{z}_t\), FHS computes risk measures by scaling the empirical quantile of those shocks by the current volatility forecast:

\[ \text{VaR}_{t+1}^{p} = \hat{\sigma}_{t+1} \cdot |\hat{z}_{(p)}| \]

Compare this with the Normal approach, which uses the same volatility but assumes the shock quantile equals \(\Phi^{-1}(p)\):

\[ \text{VaR}_{t+1}^{p,\text{Normal}} = \hat{\sigma}_{t+1} \cdot |\Phi^{-1}(p)| \]

And with standard Historical Simulation, which ignores \(\hat{\sigma}_{t+1}\) entirely and uses the percentile of raw returns:

\[ \text{VaR}_{t+1}^{p,\text{HS}} = |\text{Percentile}(\{R_{t+1-\tau}\}_{\tau=1}^{m},\; 100p)| \]

Note

How the return series is simulated. A GARCH(1,1) process is simulated with a fixed long-run daily volatility of 1.5%, using the \(\alpha\) and \(\beta\) sliders below and the chosen shock distribution. The standardized residuals are then extracted using the true GARCH parameters. The volatility forecast slider \(\hat{\sigma}_{t+1}\) is independent of the simulation: it controls only the scaling applied to the empirical shock quantile when computing VaR and ES, letting you explore how the same set of shocks translates into different risk measures at different volatility levels.

Tip

How to experiment

Switch to “Standardized t” with low degrees of freedom to see FHS capture fat tails. Then vary the volatility forecast \(\hat{\sigma}_{t+1}\): watch FHS VaR and Normal VaR scale linearly while Standard HS stays flat. In the “Exceeding history” tab, increase \(\hat{\sigma}_{t+1}\) to see FHS generate losses larger than any observed in the historical sample.

viewof fvSigma = Inputs.range([0.5, 5.0], {
  label: "σ̂(t+1) volatility forecast (%)",
  step: 0.1,
  value: 2.0
})

viewof fvConf = Inputs.range([0.5, 5.0], {
  label: "Tail probability p (%)",
  step: 0.1,
  value: 1.0
})

viewof fvAlpha = Inputs.range([0.01, 0.25], {
  label: "α (GARCH reaction)",
  step: 0.01,
  value: 0.10
})

viewof fvBeta = Inputs.range([0.50, 0.98], {
  label: "β (GARCH persistence)",
  step: 0.01,
  value: 0.85
})

viewof fvDist = Inputs.radio(["Normal", "Standardized t"], {
  label: "True shock distribution",
  value: "Normal"
})

viewof fvDf = fvDist === "Standardized t" ? Inputs.range([4, 30], {
  label: "Degrees of freedom d",
  step: 1,
  value: 8
}) : Object.assign(html`<span></span>`, {value: 8})

viewof fvN = Inputs.range([500, 3000], {
  label: "Sample size",
  step: 100,
  value: 2000
})

mutable fvSeed = 0

viewof fvReseed = {
  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 sample</button>`
  btn.onclick = () => { mutable fvSeed++ }
  return btn
}

fvVL = (1.5 / 100) ** 2  // fixed long-run vol for simulation
fvOmega = fvVL * (1 - fvAlpha - fvBeta)

fvSim = {
  fvSeed
  const rng = rng_utils.mulberry32(77 + fvSeed)
  const n = fvN
  const alpha = fvAlpha, beta = fvBeta, omega = fvOmega, vl = fvVL

  const ret = new Array(n)
  const sig2 = new Array(n)

  sig2[0] = vl
  for (let t = 0; t < n; t++) {
    const z = fvDist === "Normal"
      ? rng_utils.boxMuller(rng)
      : rng_utils.stdTSample(rng, fvDf)
    ret[t] = Math.sqrt(sig2[t]) * z
    if (t < n - 1) {
      sig2[t + 1] = omega + alpha * ret[t] * ret[t] + beta * sig2[t]
    }
  }

  // Standardized residuals
  const stdResid = ret.map((r, i) => r / Math.sqrt(sig2[i]))

  return { ret, sig2, stdResid }
}

// Compute FHS, Normal, and HS risk measures
fvResults = {
  const p = fvConf / 100
  const sigma = fvSigma / 100
  const n = fvSim.stdResid.length

  // Sort standardized shocks
  const shocks = fvSim.stdResid.slice().sort((a, b) => a - b)

  // Empirical percentile with interpolation
  const pos = fpRound(n * p)
  const lo = Math.max(0, Math.floor(pos) - 1)
  const hi = Math.min(n - 1, Math.ceil(pos) - 1)
  const frac = pos - Math.floor(pos)
  const shockQuantile = (lo === hi || frac === 0) ? shocks[lo] : shocks[lo] * (1 - frac) + shocks[hi] * frac

  // FHS VaR and ES
  const fhsVaR = sigma * Math.abs(shockQuantile)
  const tailCount = Math.max(1, Math.ceil(pos) - 1)
  const tailShocks = shocks.slice(0, tailCount)
  const shockESval = tailShocks.reduce((a, b) => a + b, 0) / tailShocks.length
  const fhsES = sigma * Math.abs(shockESval)

  // Normal VaR and ES
  const normVaR = sigma * Math.abs(qnorm(p))
  const normES = normalES(p, sigma)

  // Standard HS VaR (from raw returns)
  const rawSorted = fvSim.ret.slice().sort((a, b) => a - b)
  const rawQuantile = (lo === hi || frac === 0) ? rawSorted[lo] : rawSorted[lo] * (1 - frac) + rawSorted[hi] * frac
  const hsVaR = Math.abs(rawQuantile)
  const rawTail = rawSorted.slice(0, tailCount)
  const hsES = Math.abs(rawTail.reduce((a, b) => a + b, 0) / rawTail.length)

  // Worst loss info
  let worstIdx = 0
  for (let i = 1; i < fvSim.ret.length; i++) {
    if (fvSim.ret[i] < fvSim.ret[worstIdx]) worstIdx = i
  }
  const worstRet = fvSim.ret[worstIdx]
  const worstSigma = Math.sqrt(fvSim.sig2[worstIdx])
  const worstShock = fvSim.stdResid[worstIdx]
  const fhsImpliedLoss = sigma * Math.abs(worstShock)

  return {
    fhsVaR, fhsES, normVaR, normES, hsVaR, hsES,
    shockQuantile, shockESval, p, sigma, n, shocks,
    tailCount, worstRet, worstSigma, worstShock, worstIdx, fhsImpliedLoss
  }
}

{
  const shocks = fvResults.shocks
  const sq = fvResults.shockQuantile
  const nq = qnorm(fvResults.p)

  return Plot.plot({
    height: 380, marginLeft: 60,
    x: { label: "Standardized shock ẑ", grid: false },
    y: { label: "Count", grid: true },
    marks: [
      Plot.rectY(shocks, Plot.binX({ y: "count" }, {
        x: d => d,
        fill: d => d < sq ? "#d62728" : "#ccc",
        fillOpacity: 0.6,
        thresholds: 80
      })),
      Plot.ruleX([sq], { stroke: "#d62728", strokeWidth: 2 }),
      Plot.ruleX([nq], { stroke: "#2f71d5", strokeWidth: 2, strokeDasharray: "6 3" }),
      Plot.text([[sq, 0]], { x: d => d[0], y: d => d[1], text: d => `FHS: ${fmt(d[0], 3)}`, dy: -10, dx: -5, textAnchor: "end", fill: "#d62728", fontWeight: 700, fontSize: 12 }),
      Plot.text([[nq, 0]], { x: d => d[0], y: d => d[1], text: d => `Normal: ${fmt(d[0], 3)}`, dy: -10, dx: 5, textAnchor: "start", fill: "#2f71d5", fontWeight: 700, fontSize: 12 })
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">Histogram of standardized shocks. The <span style="color:#d62728;font-weight:700;">red line</span> is the ${fvConf}% empirical percentile used by FHS. The <span style="color:#2f71d5;font-weight:700;">blue dashed line</span> is Φ<sup>−1</sup>(${fvConf}%) = ${fmt(qnorm(fvResults.p), 3)} used by the Normal approach. The red shaded area contains the ${fvResults.tailCount} observations used to compute ES.</p>`

viewof fvVarLegend = legend([
  { key: "fhs", label: "FHS VaR", color: "#d62728", type: "line" },
  { key: "normal", label: "Normal VaR", color: "#2f71d5", type: "line" },
  { key: "hs", label: "Standard HS VaR (flat)", color: "#e67e22", type: "dashed" }
])

{
  const p = fvResults.p
  const sq = Math.abs(fvResults.shockQuantile)
  const nq = Math.abs(qnorm(p))
  const h = fvVarLegend

  // Generate sweep data
  const sweepData = []
  for (let s = 0.5; s <= 5.0; s += 0.1) {
    const sigma = s / 100
    sweepData.push({ sigma: s, FHS: sigma * sq * 100, Normal: sigma * nq * 100 })
  }

  const curSigma = fvSigma
  const hsVaR100 = fvResults.hsVaR * 100

  return Plot.plot({
    height: 380, marginLeft: 60,
    x: { label: "Volatility forecast σ̂(t+1) (%)", grid: false },
    y: { label: "VaR (% of portfolio)", grid: true },
    marks: [
      ...(!h.has("hs") ? [
        Plot.ruleY([hsVaR100], { stroke: "#e67e22", strokeWidth: 1.5, strokeDasharray: "4 3" })
      ] : []),
      ...(!h.has("fhs") ? [
        Plot.line(sweepData, { x: "sigma", y: "FHS", stroke: "#d62728", strokeWidth: 2 }),
        Plot.dot([{ x: curSigma, y: curSigma / 100 * sq * 100 }], { x: "x", y: "y", r: 5, fill: "#d62728" })
      ] : []),
      ...(!h.has("normal") ? [
        Plot.line(sweepData, { x: "sigma", y: "Normal", stroke: "#2f71d5", strokeWidth: 2 }),
        Plot.dot([{ x: curSigma, y: curSigma / 100 * nq * 100 }], { x: "x", y: "y", r: 5, fill: "#2f71d5" })
      ] : [])
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">Both FHS and Normal VaR scale linearly with σ̂<sub>t+1</sub>, but FHS has a steeper slope when shocks are fat-tailed (|ẑ<sub>(p)</sub>| > |Φ<sup>−1</sup>(p)|). Standard HS VaR is constant because it ignores the current volatility.</p>`

{
  const r = fvResults
  const worstLoss = Math.abs(r.worstRet) * 100
  const worstSig = r.worstSigma * 100
  const worstZ = Math.abs(r.worstShock)
  const fhsLoss = r.fhsImpliedLoss * 100
  const exceeds = fhsLoss > worstLoss
  const curSig = fvSigma

  return html`<div style="max-width:600px;margin:0 auto;">
    <table class="table" style="width:100%;">
    <thead><tr><th colspan="2" style="text-align:center;">Worst historical loss analysis</th></tr></thead>
    <tbody>
    <tr><td style="font-weight:500;">Worst raw loss</td>
        <td style="text-align:right;font-weight:700;">${fmt(worstLoss, 3)}%</td></tr>
    <tr><td style="font-weight:500;">Day's volatility σ<sub>${r.worstIdx + 1}</sub></td>
        <td style="text-align:right;">${fmt(worstSig, 3)}%</td></tr>
    <tr><td style="font-weight:500;">Implied shock |ẑ| = loss / σ</td>
        <td style="text-align:right;">${fmt(worstZ, 3)}</td></tr>
    <tr style="border-top:2px solid #333;"><td style="font-weight:500;">Current σ̂<sub>t+1</sub></td>
        <td style="text-align:right;">${fmt(curSig, 1)}%</td></tr>
    <tr><td style="font-weight:500;">FHS-implied loss = |ẑ| × σ̂<sub>t+1</sub></td>
        <td style="text-align:right;font-weight:700;font-size:1.1em;color:${exceeds ? '#d62728' : '#2e7d32'};">${fmt(fhsLoss, 3)}%</td></tr>
    </tbody></table>
    <div style="padding:12px;border-radius:8px;background:${exceeds ? '#fff3e0' : '#e8f5e9'};text-align:center;margin-top:8px;">
      ${exceeds
        ? html`<strong style="color:#d62728;">FHS generates a ${fmt(fhsLoss, 3)}% loss, exceeding the worst historical loss of ${fmt(worstLoss, 3)}%.</strong><br><span style="font-size:0.9rem;">A moderate past shock (${fmt(worstZ, 2)}σ) combined with high current volatility (${fmt(curSig, 1)}%) produces a loss never observed in the data.</span>`
        : html`<strong style="color:#2e7d32;">FHS-implied loss (${fmt(fhsLoss, 3)}%) is within the historical range.</strong><br><span style="font-size:0.9rem;">Current volatility is not high enough to push the worst shock beyond the observed maximum.</span>`
      }
    </div>
  </div>`
}

{
  const r = fvResults
  const fhsRatio = r.fhsES / r.fhsVaR
  const normRatio = r.normES / r.normVaR
  const hsRatio = r.hsES / r.hsVaR

  return html`<table class="table" style="width:100%;">
<thead><tr>
  <th>Measure</th>
  <th style="color:#d62728;">FHS</th>
  <th style="color:#2f71d5;">Normal</th>
  <th style="color:#e67e22;">Standard HS</th>
</tr></thead>
<tbody>
<tr><td style="font-weight:500;">VaR (% of portfolio)</td>
    <td style="font-weight:700;">${fmt(r.fhsVaR * 100, 4)}%</td>
    <td>${fmt(r.normVaR * 100, 4)}%</td>
    <td>${fmt(r.hsVaR * 100, 4)}%</td></tr>
<tr><td style="font-weight:500;">ES (% of portfolio)</td>
    <td style="font-weight:700;">${fmt(r.fhsES * 100, 4)}%</td>
    <td>${fmt(r.normES * 100, 4)}%</td>
    <td>${fmt(r.hsES * 100, 4)}%</td></tr>
<tr><td style="font-weight:500;">ES / VaR ratio</td>
    <td>${fmt(fhsRatio, 4)}</td>
    <td>${fmt(normRatio, 4)}</td>
    <td>${fmt(hsRatio, 4)}</td></tr>
<tr><td style="font-weight:500;">Shock quantile</td>
    <td>ẑ<sub>(p)</sub> = ${fmt(r.shockQuantile, 4)}</td>
    <td>Φ<sup>−1</sup>(p) = ${fmt(qnorm(r.p), 4)}</td>
    <td>n/a</td></tr>
</tbody></table>
<p style="color:#666;font-size:0.85rem;">FHS uses the empirical shock quantile ẑ<sub>(p)</sub> scaled by the volatility forecast σ̂<sub>t+1</sub> = ${fvSigma}%. The Normal approach uses Φ<sup>−1</sup>(p) instead. Standard HS uses the raw return percentile directly, ignoring current volatility. When shocks are fat-tailed, the FHS shock quantile is more negative than the normal one, producing higher VaR and ES.</p>`
}

3. Method comparison: HS vs. GARCH-Normal vs. FHS

To see why FHS matters in practice, consider a volatility regime change: a calm period followed by a stress period. How quickly does each risk measure adapt?

Standard HS uses a rolling window of raw returns. After the regime shift, it takes many days for the high-volatility observations to dominate the window, so VaR adapts slowly.
GARCH-Normal captures volatility dynamics immediately through \(\hat{\sigma}_{t+1}\), but it uses the normal quantile, which underestimates tail risk if shocks are fat-tailed.
FHS adapts immediately (through \(\hat{\sigma}_{t+1}\)) and captures fat tails (through the empirical shock quantile).

Note

How the return series is simulated. A two-regime GARCH(1,1) process is simulated with fixed \(\alpha = 0.10\) and \(\beta = 0.85\) over \(N = 2000\) days. At the regime switch point, the GARCH intercept \(\omega\) changes from \(\omega_1 = V_{L,1}(1-\alpha-\beta)\) to \(\omega_2 = V_{L,2}(1-\alpha-\beta)\), shifting the long-run variance target. The conditional variance itself transitions gradually through the GARCH recursion \(\sigma^2_{t+1} = \omega + \alpha R^2_t + \beta \sigma^2_t\) (no discontinuous jump). Shocks are always drawn from a standardized Student-t distribution with degrees of freedom \(d\). All three methods use the true GARCH conditional volatility \(\sigma_t\) (known from the simulation), so the comparison isolates the effect of the shock distribution assumption (normal quantile vs. empirical quantile) and the volatility responsiveness (rolling window vs. GARCH forecast).

Note

In-sample illustration. This comparison is performed in-sample: the same data used to build the shock database is also used to evaluate exceedances. In practice, model performance should be assessed out-of-sample, applying VaR forecasts to future returns not used in estimation. In-sample results tend to be optimistic, particularly for FHS, which benefits from seeing the shocks it is evaluated against. The purpose here is to illustrate the qualitative differences between methods, not to provide a realistic backtesting exercise.

Tip

How to experiment

Set a large contrast between calm and stress volatilities (e.g., 1% vs. 4%). After the regime switch, watch how HS VaR lags behind while FHS and GARCH-Normal adapt quickly. Then check the “Exceedance rates” tab: HS will be too conservative before the switch and too lenient after, while FHS should maintain rates closer to the target \(p\). Lower the degrees of freedom to increase the gap between GARCH-Normal and FHS.

viewof mcCalmVol = Inputs.range([0.5, 2.0], {
  label: "Calm-period daily volatility (%)",
  step: 0.1,
  value: 1.0
})

viewof mcStressVol = Inputs.range([2.0, 5.0], {
  label: "Stress-period daily volatility (%)",
  step: 0.1,
  value: 3.0
})

viewof mcSwitch = Inputs.range([20, 80], {
  label: "Regime switch point (%)",
  step: 5,
  value: 60
})

viewof mcConf = Inputs.range([0.5, 5.0], {
  label: "Tail probability p (%)",
  step: 0.1,
  value: 1.0
})

viewof mcDf = Inputs.range([4, 30], {
  label: "Degrees of freedom d",
  step: 1,
  value: 5
})

viewof mcWindow = Inputs.range([100, 500], {
  label: "HS rolling window",
  step: 50,
  value: 250
})

mutable mcSeed = 0

viewof mcReseed = {
  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 sample</button>`
  btn.onclick = () => { mutable mcSeed++ }
  return btn
}

// Two-regime GARCH(1,1) simulation
mcN = 2000

mcSim = {
  mcSeed
  const rng = rng_utils.mulberry32(42 + mcSeed)
  const N = mcN
  const switchIdx = Math.floor(N * mcSwitch / 100)
  const alpha = 0.10, beta = 0.85

  const vl1 = (mcCalmVol / 100) ** 2
  const vl2 = (mcStressVol / 100) ** 2
  const omega1 = vl1 * (1 - alpha - beta)
  const omega2 = vl2 * (1 - alpha - beta)

  const ret = new Array(N)
  const sig2 = new Array(N)
  const shock = new Array(N)

  sig2[0] = vl1
  for (let t = 0; t < N; t++) {
    const z = rng_utils.stdTSample(rng, mcDf)
    shock[t] = z
    ret[t] = Math.sqrt(sig2[t]) * z
    if (t < N - 1) {
      const omega = (t + 1 < switchIdx) ? omega1 : omega2
      sig2[t + 1] = omega + alpha * ret[t] * ret[t] + beta * sig2[t]
    }
  }

  return { ret, sig2, shock, switchIdx, alpha, beta }
}

// Rolling VaR computation for all three methods
mcRolling = {
  const p = mcConf / 100
  const { ret, sig2, shock, switchIdx } = mcSim
  const m = mcWindow
  const N = mcN
  const results = []

  for (let t = m; t < N; t++) {
    const sigT = Math.sqrt(sig2[t])

    // Standard HS: percentile of raw returns in rolling window
    const rawWindow = ret.slice(t - m, t).slice().sort((a, b) => a - b)
    const pos = fpRound(m * p)
    const lo = Math.max(0, Math.floor(pos) - 1)
    const hi = Math.min(m - 1, Math.ceil(pos) - 1)
    const frac = pos - Math.floor(pos)
    const rawQ = (lo === hi || frac === 0) ? rawWindow[lo] : rawWindow[lo] * (1 - frac) + rawWindow[hi] * frac
    const hsVaR = Math.abs(rawQ)

    // GARCH-Normal
    const gnVaR = sigT * Math.abs(qnorm(p))

    // FHS: expanding window of all shocks up to t
    const shockWindow = shock.slice(0, t).slice().sort((a, b) => a - b)
    const nShocks = shockWindow.length
    const sPos = fpRound(nShocks * p)
    const sLo = Math.max(0, Math.floor(sPos) - 1)
    const sHi = Math.min(nShocks - 1, Math.ceil(sPos) - 1)
    const sFrac = sPos - Math.floor(sPos)
    const shockQ = (sLo === sHi || sFrac === 0) ? shockWindow[sLo] : shockWindow[sLo] * (1 - sFrac) + shockWindow[sHi] * sFrac
    const fhsVaR = sigT * Math.abs(shockQ)

    // Exceedance check
    const loss = -ret[t]
    results.push({
      t: t + 1,
      ret: ret[t],
      loss,
      hsVaR, gnVaR, fhsVaR,
      hsExceed: loss > hsVaR,
      gnExceed: loss > gnVaR,
      fhsExceed: loss > fhsVaR,
      preSwitch: t < switchIdx
    })
  }
  return results
}

// Exceedance statistics
mcStats = {
  const pre = mcRolling.filter(d => d.preSwitch)
  const post = mcRolling.filter(d => !d.preSwitch)
  const all = mcRolling
  const p = mcConf / 100

  const count = (arr, key) => arr.filter(d => d[key]).length
  const rate = (arr, key) => arr.length > 0 ? count(arr, key) / arr.length : 0

  return {
    p,
    total: all.length,
    preN: pre.length,
    postN: post.length,
    hs: {
      pre: rate(pre, "hsExceed"), post: rate(post, "hsExceed"), all: rate(all, "hsExceed"),
      preCount: count(pre, "hsExceed"), postCount: count(post, "hsExceed"), allCount: count(all, "hsExceed")
    },
    gn: {
      pre: rate(pre, "gnExceed"), post: rate(post, "gnExceed"), all: rate(all, "gnExceed"),
      preCount: count(pre, "gnExceed"), postCount: count(post, "gnExceed"), allCount: count(all, "gnExceed")
    },
    fhs: {
      pre: rate(pre, "fhsExceed"), post: rate(post, "fhsExceed"), all: rate(all, "fhsExceed"),
      preCount: count(pre, "fhsExceed"), postCount: count(post, "fhsExceed"), allCount: count(all, "fhsExceed")
    }
  }
}

viewof mcVarLegend = legend([
  { key: "returns", label: "Returns", color: "#bbb", type: "line" },
  { key: "hs", label: "Standard HS VaR", color: "#e67e22", type: "line" },
  { key: "gn", label: "GARCH-Normal VaR", color: "#2f71d5", type: "line" },
  { key: "fhs", label: "FHS VaR", color: "#d62728", type: "line" }
])

{
  const switchDay = mcSim.switchIdx + 1
  const h = mcVarLegend

  return Plot.plot({
    height: 420, marginLeft: 60,
    x: { label: "Day", grid: false },
    y: { label: "Return / VaR", grid: true, tickFormat: pctFmt },
    marks: [
      Plot.ruleY([0], { stroke: "#eee" }),
      Plot.ruleX([switchDay], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
      ...(!h.has("returns") ? [
        Plot.line(mcRolling, { x: "t", y: "ret", stroke: "#bbb", strokeWidth: 0.4 })
      ] : []),
      ...(!h.has("hs") ? [
        Plot.line(mcRolling, { x: "t", y: d => -d.hsVaR, stroke: "#e67e22", strokeWidth: 1.5 })
      ] : []),
      ...(!h.has("gn") ? [
        Plot.line(mcRolling, { x: "t", y: d => -d.gnVaR, stroke: "#2f71d5", strokeWidth: 1.5 })
      ] : []),
      ...(!h.has("fhs") ? [
        Plot.line(mcRolling, { x: "t", y: d => -d.fhsVaR, stroke: "#d62728", strokeWidth: 1.5 })
      ] : []),
      Plot.text([[switchDay, 0]], { x: d => d[0], y: d => d[1], text: ["Regime switch"], dy: -5, textAnchor: "middle", fontSize: 11, fill: "#666" })
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">VaR lines are plotted as negative values (loss thresholds). After the regime switch, Standard HS (orange) adjusts slowly because it relies on a rolling window. GARCH-Normal (blue) and FHS (red) adapt immediately. FHS exceeds GARCH-Normal when shocks are fat-tailed.</p>`

viewof mcExcLegend = legend([
  { key: "returns", label: "Returns", color: "#ddd", type: "line" },
  { key: "hs", label: "HS exceedance", color: "#e67e22", type: "dot" },
  { key: "gn", label: "GARCH-Normal exceedance", color: "#2f71d5", type: "dot" },
  { key: "fhs", label: "FHS exceedance", color: "#d62728", type: "dot" }
])

{
  const switchDay = mcSim.switchIdx + 1
  const h = mcExcLegend

  // Build exceedance data with method labels and keys
  const excData = []
  for (const d of mcRolling) {
    if (d.hsExceed) excData.push({ t: d.t, ret: d.ret, method: "HS", key: "hs", color: "#e67e22" })
    if (d.gnExceed) excData.push({ t: d.t, ret: d.ret, method: "GARCH-Normal", key: "gn", color: "#2f71d5" })
    if (d.fhsExceed) excData.push({ t: d.t, ret: d.ret, method: "FHS", key: "fhs", color: "#d62728" })
  }
  const visibleExc = excData.filter(d => !h.has(d.key))

  return Plot.plot({
    height: 420, marginLeft: 60,
    x: { label: "Day", grid: false },
    y: { label: "Return", grid: true, tickFormat: pctFmt },
    marks: [
      Plot.ruleY([0], { stroke: "#eee" }),
      Plot.ruleX([switchDay], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
      ...(!h.has("returns") ? [
        Plot.line(mcRolling, { x: "t", y: "ret", stroke: "#ddd", strokeWidth: 0.4 })
      ] : []),
      ...(visibleExc.length > 0 ? [
        Plot.dot(visibleExc, { x: "t", y: "ret", fill: "color", r: 3, fillOpacity: 0.7 })
      ] : [])
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">Each dot marks a day where the actual loss exceeded the method's VaR. A cluster of exceedances after the regime switch indicates the method was too slow to react. Fewer exceedances overall (closer to p × N) indicates better calibration.</p>`

{
  const s = mcStats
  const fmtPct = x => (x * 100).toFixed(2) + "%"
  const highlight = (rate, target) => {
    const ratio = rate / target
    if (ratio > 2) return 'color:#d62728;font-weight:700;'
    if (ratio < 0.5) return 'color:#2e7d32;font-weight:700;'
    return ''
  }

  return html`<table class="table" style="width:100%;">
<thead><tr>
  <th>Method</th>
  <th>Pre-switch rate</th>
  <th>Post-switch rate</th>
  <th>Overall rate</th>
  <th>Target p</th>
</tr></thead>
<tbody>
<tr><td style="font-weight:500;color:#e67e22;">Standard HS</td>
    <td style="${highlight(s.hs.pre, s.p)}">${fmtPct(s.hs.pre)} (${s.hs.preCount}/${s.preN})</td>
    <td style="${highlight(s.hs.post, s.p)}">${fmtPct(s.hs.post)} (${s.hs.postCount}/${s.postN})</td>
    <td style="${highlight(s.hs.all, s.p)}">${fmtPct(s.hs.all)} (${s.hs.allCount}/${s.total})</td>
    <td>${fmtPct(s.p)}</td></tr>
<tr><td style="font-weight:500;color:#2f71d5;">GARCH-Normal</td>
    <td style="${highlight(s.gn.pre, s.p)}">${fmtPct(s.gn.pre)} (${s.gn.preCount}/${s.preN})</td>
    <td style="${highlight(s.gn.post, s.p)}">${fmtPct(s.gn.post)} (${s.gn.postCount}/${s.postN})</td>
    <td style="${highlight(s.gn.all, s.p)}">${fmtPct(s.gn.all)} (${s.gn.allCount}/${s.total})</td>
    <td>${fmtPct(s.p)}</td></tr>
<tr><td style="font-weight:500;color:#d62728;">FHS</td>
    <td style="${highlight(s.fhs.pre, s.p)}">${fmtPct(s.fhs.pre)} (${s.fhs.preCount}/${s.preN})</td>
    <td style="${highlight(s.fhs.post, s.p)}">${fmtPct(s.fhs.post)} (${s.fhs.postCount}/${s.postN})</td>
    <td style="${highlight(s.fhs.all, s.p)}">${fmtPct(s.fhs.all)} (${s.fhs.allCount}/${s.total})</td>
    <td>${fmtPct(s.p)}</td></tr>
</tbody></table>
<p style="color:#666;font-size:0.85rem;">Rates highlighted in <span style="color:#d62728;font-weight:700;">red</span> exceed twice the target (too many exceedances, VaR too low). Rates in <span style="color:#2e7d32;font-weight:700;">green</span> are below half the target (too few exceedances, VaR too conservative). A well-calibrated model should have rates close to p = ${fmtPct(s.p)} in both periods.</p>`
}

{
  const s = mcStats
  const expected = s.total * s.p

  const barData = [
    { method: "Standard HS", count: s.hs.allCount },
    { method: "GARCH-Normal", count: s.gn.allCount },
    { method: "FHS", count: s.fhs.allCount }
  ]

  return Plot.plot({
    height: 350, marginLeft: 120, marginRight: 40,
    x: { label: "Exceedance count", grid: true },
    y: { label: null, domain: ["Standard HS", "GARCH-Normal", "FHS"] },
    color: { domain: ["Standard HS", "GARCH-Normal", "FHS"], range: ["#e67e22", "#2f71d5", "#d62728"] },
    marks: [
      Plot.barX(barData, { x: "count", y: "method", fill: "method" }),
      Plot.ruleX([expected], { stroke: "#333", strokeWidth: 2, strokeDasharray: "6 3" }),
      Plot.text(barData, { x: "count", y: "method", text: d => d.count.toString(), dx: 8, textAnchor: "start", fontWeight: 700, fontSize: 13 }),
      Plot.text([[expected, "FHS"]], { x: d => d[0], y: d => d[1], text: [`Expected: ${fmt(expected, 1)}`], dy: 40, textAnchor: "middle", fontSize: 11, fill: "#333" })
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">Total exceedance counts for each method over ${mcStats.total} evaluation days. The dashed line marks the expected count under perfect calibration (p × N = ${fmt(mcStats.total * mcStats.p, 1)}). A well-calibrated method should have a count close to this line.</p>`

References

Barone-Adesi, Giovanni, Kostas Giannopoulos, and Les Vosper. 1999. “VaR Without Correlations for Portfolios of Derivative Securities.” Journal of Futures Markets 19 (5): 583–602.

Christoffersen, Peter F. 2012. Elements of Financial Risk Management. 2nd ed. Academic Press.