GARCH Volatility Models

Interactive exploration of EWMA and GARCH models for time-varying conditional volatility

Financial return volatility is not constant: it clusters in time and mean-reverts to a long-run level. This page explores the two workhorse models for tracking conditional volatility: the Exponentially Weighted Moving Average (EWMA/RiskMetrics) and the Generalized Autoregressive Conditional Heteroskedasticity (GARCH, Bollerslev 1986) family. Both models decompose daily returns as

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

where \(\sigma_{t+1}\) is the conditional standard deviation forecast and \(z_{t+1}\) is the standardized innovation. The goal is to model and forecast \(\sigma^2_{t+1}\) (see Hull 2023, secs. 8.6–8.10; Christoffersen 2012, chap. 4).

// ============================================================
// SHARED UTILITIES
// ============================================================

// Seeded PRNG (Mulberry32)
prng = {
  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)
  }
  return { mulberry32, boxMuller }
}

// 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) + "%"

1. EWMA vs GARCH(1,1)

The EWMA model updates variance as a weighted average of yesterday’s variance and squared return, with all weight on recent data and no mean reversion:

\[ \sigma^2_{t+1} = \lambda\,\sigma^2_t + (1-\lambda)\,R^2_t, \quad \lambda = 0.94 \]

The GARCH(1,1) model adds a pull toward long-run variance \(V_L = \omega/(1-\alpha-\beta)\):

\[ \sigma^2_{t+1} = \omega + \alpha\,R^2_t + \beta\,\sigma^2_t, \quad \alpha + \beta < 1 \]

Both models are applied to the same simulated return series (generated from the GARCH process) so that the comparison is like-for-like.

Note

What happens when \(\alpha + \beta \geq 1\)?

The stationarity condition \(\alpha + \beta < 1\) ensures that the long-run variance \(V_L = \omega/(1-\alpha-\beta)\) is well-defined and positive, and that shocks to variance eventually decay. When this condition is violated:

\(\alpha + \beta = 1\) (the EWMA/RiskMetrics case): \(\omega = 0\) and the long-run variance is undefined. Variance follows a random walk: any shock persists forever in the forecast, and the model has no “anchor” to revert to. This is the EWMA model with \(\lambda = \beta\).
\(\alpha + \beta > 1\): The process is explosive. Shocks to variance are amplified over time rather than decaying, and the unconditional variance is infinite. The model is mean-averting rather than mean-reverting.

In practice, estimated GARCH models for financial data almost always yield \(\alpha + \beta\) slightly below 1 (typically 0.95 to 0.99), confirming that variance is highly persistent but ultimately stationary.

Tip

How to experiment

Adjust the GARCH parameters \(\alpha\) and \(\beta\) to control persistence and reactivity. Then compare with EWMA by varying \(\lambda\). After a volatility spike, watch how GARCH variance reverts to the long-run level (dashed line) while EWMA stays elevated. Check the ACF tabs to see which model better removes autocorrelation from squared returns. Try pushing \(\alpha + \beta\) close to or above 1 to see the stationarity constraint in action.

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

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

viewof cmpLambda = Inputs.range([0.88, 0.99], {
  label: "λ (EWMA decay factor)",
  step: 0.01,
  value: 0.94
})

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

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

mutable cmpSeed = 0

viewof cmpReseed = {
  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 cmpSeed++ }
  return btn
}

cmpPersistence = cmpAlpha + cmpBeta
cmpVL = (cmpLRVol / 100) ** 2
cmpOmega = cmpVL * (1 - cmpAlpha - cmpBeta)
cmpStationary = cmpPersistence < 1

// Simulate returns from GARCH(1,1) and compute both GARCH and EWMA variance paths
cmpSim = {
  cmpSeed  // reactivity trigger
  const rng = prng.mulberry32(42 + cmpSeed)
  const n = cmpN
  const alpha = cmpAlpha, beta = cmpBeta, omega = cmpOmega, vl = cmpVL, lambda = cmpLambda

  // Arrays
  const ret = new Array(n)
  const gVar = new Array(n)   // GARCH conditional variance
  const eVar = new Array(n)   // EWMA conditional variance

  // Initialize
  gVar[0] = vl
  eVar[0] = vl

  for (let t = 0; t < n; t++) {
    const z = prng.boxMuller(rng)
    ret[t] = Math.sqrt(gVar[t]) * z

    if (t < n - 1) {
      gVar[t + 1] = omega + alpha * ret[t] * ret[t] + beta * gVar[t]
      eVar[t + 1] = lambda * eVar[t] + (1 - lambda) * ret[t] * ret[t]
    }
  }

  return { ret, gVar, eVar }
}

// Prepare time series data for plotting
cmpPlotData = {
  const { ret, gVar, eVar } = cmpSim
  return ret.map((r, i) => ({
    t: i + 1,
    ret: r,
    retSq: r * r,
    gVol: Math.sqrt(gVar[i]),
    eVol: Math.sqrt(eVar[i]),
    gVar: gVar[i],
    eVar: eVar[i],
    gStdSq: r * r / gVar[i],
    eStdSq: r * r / eVar[i]
  }))
}

cmpACFLags = 50
cmpACFRetSq = computeACF(cmpSim.ret.map(r => r * r), cmpACFLags)
cmpACFGStdSq = computeACF(cmpPlotData.map(d => d.gStdSq), cmpACFLags)
cmpACFEStdSq = computeACF(cmpPlotData.map(d => d.eStdSq), cmpACFLags)
cmpBartlett = 2 / Math.sqrt(cmpN)

cmpLBRetSq = ljungBox(cmpSim.ret.map(r => r * r), 15)
cmpLBGStdSq = ljungBox(cmpPlotData.map(d => d.gStdSq), 15)
cmpLBEStdSq = ljungBox(cmpPlotData.map(d => d.eStdSq), 15)

Plot.plot({
  height: 350,
  marginLeft: 60,
  x: { label: "Day", grid: false },
  y: { label: "Return", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.line(cmpPlotData, { x: "t", y: "ret", stroke: "#333", strokeWidth: 0.5 })
  ]
})

html`<p style="color:#666;font-size:0.85rem;">Simulated daily returns from GARCH(1,1) with α = ${cmpAlpha}, β = ${cmpBeta}, persistence = ${fmt(cmpPersistence, 2)}. Long-run daily volatility = ${cmpLRVol}%.</p>`

Plot.plot({
  height: 380,
  marginLeft: 60,
  x: { label: "Day", grid: false },
  y: { label: "Daily volatility (σ)", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([cmpLRVol / 100], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
    Plot.line(cmpPlotData, { x: "t", y: "eVol", stroke: "#d62728", strokeWidth: 1.5 }),
    Plot.line(cmpPlotData, { x: "t", y: "gVol", stroke: "#2f71d5", strokeWidth: 1.5 })
  ]
})

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#2f71d5" stroke-width="2"/></svg> GARCH(1,1)</span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#d62728" stroke-width="2"/></svg> EWMA (λ = ${cmpLambda})</span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#999" stroke-width="1.5" stroke-dasharray="6 3"/></svg> Long-run volatility (${cmpLRVol}%)</span>
</div>`

html`<p style="color:#666;font-size:0.85rem;">Both models track volatility similarly in the short run. After a spike, the GARCH variance reverts toward the long-run level while EWMA stays elevated. ${!cmpStationary ? '<strong style="color:#d62728;">Warning: α + β ≥ 1, GARCH is non-stationary (equivalent to EWMA).</strong>' : ''}</p>`

Plot.plot({
  height: 350,
  marginLeft: 60,
  x: { label: "Lag", grid: false, ticks: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50] },
  y: { label: "Autocorrelation", grid: true, domain: [-0.15, Math.max(0.5, ...cmpACFRetSq.map(d => d.acf) ) * 1.1] },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.ruleY([cmpBartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.ruleY([-cmpBartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.rectY(cmpACFRetSq, { x1: d => d.lag - 0.4, x2: d => d.lag + 0.4, y: "acf", fill: "#4682b4" })
  ]
})

html`<p style="color:#666;font-size:0.85rem;">Autocorrelation of squared returns R²<sub>t</sub>. Positive autocorrelations at short lags are a signature of <strong>volatility clustering</strong>. The dashed lines show 95% Bartlett confidence bands (±${fmt(cmpBartlett, 3)}).</p>`

{
  // Offset GARCH and EWMA bars slightly for side-by-side display
  const garchBars = cmpACFGStdSq.map(d => ({ lag: d.lag, acf: d.acf }))
  const ewmaBars = cmpACFEStdSq.map(d => ({ lag: d.lag, acf: d.acf }))

  return Plot.plot({
    height: 350,
    marginLeft: 60,
    x: {
      label: "Lag",
      grid: false,
      ticks: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50]
    },
    y: { label: "Autocorrelation", grid: true, domain: [-0.15, 0.15] },
    marks: [
      Plot.ruleY([0], { stroke: "#ccc" }),
      Plot.ruleY([cmpBartlett], { stroke: "#999", strokeDasharray: "4 3" }),
      Plot.ruleY([-cmpBartlett], { stroke: "#999", strokeDasharray: "4 3" }),
      Plot.rectY(garchBars, { x1: d => d.lag - 0.4, x2: d => d.lag, y: "acf", fill: "#2f71d5" }),
      Plot.rectY(ewmaBars, { x1: d => d.lag, x2: d => d.lag + 0.4, y: "acf", fill: "#d62728" })
    ]
  })
}

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="12" height="12"><rect width="12" height="12" fill="#2f71d5"/></svg> GARCH</span>
  <span><svg width="12" height="12"><rect width="12" height="12" fill="#d62728"/></svg> EWMA</span>
</div>
<p style="color:#666;font-size:0.85rem;">Autocorrelation of standardized squared returns R²<sub>t</sub> / σ²<sub>t</sub>. A well-specified model should remove the autocorrelation, leaving values within the 95% bands.</p>`

{
  const retStats = descStats(cmpSim.ret)
  const gResid = cmpSim.ret.map((r, i) => r / Math.sqrt(cmpSim.gVar[i]))
  const eResid = cmpSim.ret.map((r, i) => r / Math.sqrt(cmpSim.eVar[i]))
  const gStats = descStats(gResid)
  const eStats = descStats(eResid)

  // Theoretical unconditional kurtosis for GARCH(1,1)
  const ab = cmpAlpha + cmpBeta
  const denom = 1 - 2 * cmpAlpha * cmpAlpha - ab * ab
  const thKurt = denom > 0 ? 3 * (1 - ab * ab) / denom : Infinity

  return html`<table style="font-size:0.9rem;border-collapse:collapse;width:100%;">
    <thead>
      <tr style="border-bottom:2px solid #333;">
        <th style="text-align:left;padding:4px 8px;">Statistic</th>
        <th style="text-align:right;padding:4px 8px;">Raw returns</th>
        <th style="text-align:right;padding:4px 8px;">GARCH std. residuals</th>
        <th style="text-align:right;padding:4px 8px;">EWMA std. residuals</th>
      </tr>
    </thead>
    <tbody>
      <tr><td style="padding:4px 8px;">Mean</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(retStats.mean, 6)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(gStats.mean, 4)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(eStats.mean, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Std. deviation</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(retStats.sd, 6)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(gStats.sd, 4)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(eStats.sd, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Skewness</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(retStats.skew, 4)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(gStats.skew, 4)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(eStats.skew, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Kurtosis</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(retStats.kurt, 4)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(gStats.kurt, 4)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(eStats.kurt, 4)}</td></tr>
      <tr style="border-top:1px solid #ccc;">
        <td style="padding:4px 8px;">Ljung-Box(15) of R²<sub>t</sub> / σ²<sub>t</sub></td>
        <td style="text-align:right;padding:4px 8px;">${fmt(cmpLBRetSq, 2)} (raw R²<sub>t</sub>)</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(cmpLBGStdSq, 2)}</td>
        <td style="text-align:right;padding:4px 8px;">${fmt(cmpLBEStdSq, 2)}</td></tr>
    </tbody>
    <tfoot>
      <tr style="border-top:2px solid #333;">
        <td colspan="4" style="padding:4px 8px;font-size:0.82rem;color:#666;">
          GARCH persistence = ${fmt(cmpPersistence, 2)}.
          Long-run variance = ${fmt(cmpVL * 10000, 4)} (×10⁻⁴).
          Theoretical kurtosis = ${thKurt === Infinity ? "∞" : fmt(thKurt, 2)}.
          Ljung-Box critical value at 95% for K=15 is 25.0.
        </td>
      </tr>
    </tfoot>
  </table>`
}

2. GARCH(1,1) simulator

The GARCH(1,1) model has three parameters (\(\omega\), \(\alpha\), \(\beta\)), but these map to interpretable quantities:

\[ \sigma^2_{t+1} = \underbrace{(1-\alpha-\beta)}_{\gamma}\,V_L + \alpha\,R^2_t + \beta\,\sigma^2_t \]

where \(V_L = \omega/(1-\alpha-\beta)\) is the long-run variance and \(\alpha + \beta\) is the persistence measuring how long a volatility shock takes to decay.

Tip

How to experiment

Increase \(\alpha\) to make volatility more reactive to recent shocks. Increase \(\beta\) to make it more persistent. Watch the conditional volatility tab and ACF diagnostics. With very high persistence (\(\alpha + \beta > 0.99\)), the model behaves almost like EWMA.

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

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

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

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

mutable g11Seed = 0

viewof g11Reseed = {
  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 g11Seed++ }
  return btn
}

g11Pers = g11Alpha + g11Beta
g11VL = (g11LRVol / 100) ** 2
g11Omega = g11VL * (1 - g11Alpha - g11Beta)
g11Stationary = g11Pers < 1

g11Sim = {
  g11Seed
  const rng = prng.mulberry32(123 + g11Seed)
  const n = g11N
  const alpha = g11Alpha, beta = g11Beta, omega = g11Omega, vl = g11VL

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

  gVar[0] = vl
  for (let t = 0; t < n; t++) {
    const z = prng.boxMuller(rng)
    ret[t] = Math.sqrt(gVar[t]) * z
    if (t < n - 1) gVar[t + 1] = omega + alpha * ret[t] * ret[t] + beta * gVar[t]
  }
  return { ret, gVar }
}

g11Data = g11Sim.ret.map((r, i) => ({
  t: i + 1, ret: r, retSq: r * r,
  gVol: Math.sqrt(g11Sim.gVar[i]),
  stdSq: r * r / g11Sim.gVar[i]
}))

g11ACFRetSq = computeACF(g11Sim.ret.map(r => r * r), 50)
g11ACFStdSq = computeACF(g11Data.map(d => d.stdSq), 50)
g11Bartlett = 2 / Math.sqrt(g11N)

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Day" }, y: { label: "Return", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.line(g11Data, { x: "t", y: "ret", stroke: "#333", strokeWidth: 0.5 })
  ]
})

html`<p style="color:#666;font-size:0.85rem;">GARCH(1,1) simulated returns. α = ${g11Alpha}, β = ${g11Beta}, persistence = ${fmt(g11Pers, 2)}, ω = ${fmt(g11Omega, 8)}.</p>`

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Day" }, y: { label: "Daily volatility (σ)", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([g11LRVol / 100], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
    Plot.line(g11Data, { x: "t", y: "gVol", stroke: "#2f71d5", strokeWidth: 1 })
  ]
})

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#2f71d5" stroke-width="2"/></svg> GARCH(1,1) σ<sub>t</sub></span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#999" stroke-width="1.5" stroke-dasharray="6 3"/></svg> Long-run volatility (${g11LRVol}%)</span>
</div>`

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Lag", ticks: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50] },
  y: { label: "Autocorrelation", grid: true, domain: [-0.15, Math.max(0.5, ...g11ACFRetSq.map(d => d.acf)) * 1.1] },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.ruleY([g11Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.ruleY([-g11Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.rectY(g11ACFRetSq, { x1: d => d.lag - 0.4, x2: d => d.lag + 0.4, y: "acf", fill: "#4682b4" })
  ]
})

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Lag", ticks: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50] },
  y: { label: "Autocorrelation", grid: true, domain: [-0.15, 0.15] },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.ruleY([g11Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.ruleY([-g11Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.rectY(g11ACFStdSq, { x1: d => d.lag - 0.4, x2: d => d.lag + 0.4, y: "acf", fill: "#2f71d5" })
  ]
})

html`<p style="color:#666;font-size:0.85rem;">The GARCH model should remove autocorrelation from standardized squared returns. Values within the dashed 95% bands indicate adequate fit.</p>`

{
  const retStats = descStats(g11Sim.ret)
  const resid = g11Sim.ret.map((r, i) => r / Math.sqrt(g11Sim.gVar[i]))
  const rStats = descStats(resid)
  const ab = g11Alpha + g11Beta
  const denom = 1 - 2 * g11Alpha * g11Alpha - ab * ab
  const thKurt = denom > 0 ? 3 * (1 - ab * ab) / denom : Infinity
  const lb15 = ljungBox(g11Data.map(d => d.stdSq), 15)

  return html`<table style="font-size:0.9rem;border-collapse:collapse;width:100%;">
    <thead>
      <tr style="border-bottom:2px solid #333;">
        <th style="text-align:left;padding:4px 8px;">Statistic</th>
        <th style="text-align:right;padding:4px 8px;">Returns</th>
        <th style="text-align:right;padding:4px 8px;">Std. residuals (z<sub>t</sub>)</th>
      </tr>
    </thead>
    <tbody>
      <tr><td style="padding:4px 8px;">Mean</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.mean, 6)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.mean, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Std. deviation</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.sd, 6)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.sd, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Skewness</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.skew, 4)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.skew, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Kurtosis</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.kurt, 2)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.kurt, 4)}</td></tr>
      <tr style="border-top:1px solid #ccc;">
        <td style="padding:4px 8px;">Ljung-Box(15) of R²/σ²</td><td style="text-align:right;padding:4px 8px;"></td><td style="text-align:right;padding:4px 8px;">${fmt(lb15, 2)}</td></tr>
    </tbody>
    <tfoot>
      <tr style="border-top:2px solid #333;">
        <td colspan="3" style="padding:4px 8px;font-size:0.82rem;color:#666;">
          Persistence = ${fmt(ab, 2)}. Long-run vol = ${g11LRVol}%/day.
          Theoretical unconditional kurtosis = ${thKurt === Infinity ? "∞ (4th moment does not exist)" : fmt(thKurt, 2)}.
          LB critical value (95%, K=15) = 25.0.
        </td>
      </tr>
    </tfoot>
  </table>`
}

3. GARCH(2,2) simulator

The GARCH(2,2) model extends GARCH(1,1) with two lags each of squared returns and variance, allowing richer autocorrelation dynamics:

\[ \sigma^2_{t+1} = \omega + \alpha_1 R^2_t + \alpha_2 R^2_{t-1} + \beta_1 \sigma^2_t + \beta_2 \sigma^2_{t-1} \]

The long-run variance is \(V_L = \omega / (1 - \alpha_1 - \alpha_2 - \beta_1 - \beta_2)\). This model can equivalently be written as a component GARCH with interpretable short-run and long-run variance factors (see Christoffersen 2012, sec. 5.3 and appendix A).

Tip

How to experiment

Compare with GARCH(1,1) using similar total persistence (\(\alpha_1 + \alpha_2 + \beta_1 + \beta_2 \approx \alpha + \beta\)). GARCH(2,2) can capture slower ACF decay in squared returns. Watch for non-stationarity when the sum of all parameters approaches 1.

viewof g22A1 = Inputs.range([0.00, 0.20], { label: "α₁", step: 0.01, value: 0.05 })

viewof g22A2 = Inputs.range([0.00, 0.15], { label: "α₂", step: 0.01, value: 0.05 })

viewof g22B1 = Inputs.range([0.10, 0.90], { label: "β₁", step: 0.01, value: 0.50 })

viewof g22B2 = Inputs.range([0.00, 0.50], { label: "β₂", step: 0.01, value: 0.30 })

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

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

mutable g22Seed = 0

viewof g22Reseed = {
  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 g22Seed++ }
  return btn
}

g22Pers = g22A1 + g22A2 + g22B1 + g22B2
g22VL = (g22LRVol / 100) ** 2
g22Omega = g22VL * Math.max(0.001, 1 - g22Pers)
g22Stationary = g22Pers < 1

g22Sim = {
  g22Seed
  const rng = prng.mulberry32(777 + g22Seed)
  const n = g22N
  const a1 = g22A1, a2 = g22A2, b1 = g22B1, b2 = g22B2
  const omega = g22Omega, vl = g22VL

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

  // Initialize first two periods
  gVar[0] = vl; gVar[1] = vl
  ret[0] = Math.sqrt(gVar[0]) * prng.boxMuller(rng)

  for (let t = 1; t < n; t++) {
    const z = prng.boxMuller(rng)
    ret[t] = Math.sqrt(gVar[t]) * z
    if (t < n - 1) {
      const r2_t = ret[t] * ret[t]
      const r2_tm1 = ret[t - 1] * ret[t - 1]
      gVar[t + 1] = omega + a1 * r2_t + a2 * r2_tm1 + b1 * gVar[t] + b2 * gVar[Math.max(0, t - 1)]
      // Floor to prevent negative variance from numerical issues
      if (gVar[t + 1] < 1e-12) gVar[t + 1] = 1e-12
    }
  }
  return { ret, gVar }
}

g22Data = g22Sim.ret.map((r, i) => ({
  t: i + 1, ret: r, retSq: r * r,
  gVol: Math.sqrt(g22Sim.gVar[i]),
  stdSq: r * r / g22Sim.gVar[i]
}))

g22ACFRetSq = computeACF(g22Sim.ret.map(r => r * r), 50)
g22ACFStdSq = computeACF(g22Data.map(d => d.stdSq), 50)
g22Bartlett = 2 / Math.sqrt(g22N)

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Day" }, y: { label: "Return", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.line(g22Data, { x: "t", y: "ret", stroke: "#333", strokeWidth: 0.5 })
  ]
})

html`<p style="color:#666;font-size:0.85rem;">GARCH(2,2) simulated returns. α₁ = ${g22A1}, α₂ = ${g22A2}, β₁ = ${g22B1}, β₂ = ${g22B2}. Persistence = ${fmt(g22Pers, 2)}. ${!g22Stationary ? '<strong style="color:#d62728;">Non-stationary!</strong>' : ''}</p>`

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Day" }, y: { label: "Daily volatility (σ)", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([g22LRVol / 100], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
    Plot.line(g22Data, { x: "t", y: "gVol", stroke: "#e67e22", strokeWidth: 1 })
  ]
})

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#e67e22" stroke-width="2"/></svg> GARCH(2,2) σ<sub>t</sub></span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#999" stroke-width="1.5" stroke-dasharray="6 3"/></svg> Long-run volatility (${g22LRVol}%)</span>
</div>`

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Lag", ticks: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50] },
  y: { label: "Autocorrelation", grid: true, domain: [-0.15, Math.max(0.5, ...g22ACFRetSq.map(d => d.acf)) * 1.1] },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.ruleY([g22Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.ruleY([-g22Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.rectY(g22ACFRetSq, { x1: d => d.lag - 0.4, x2: d => d.lag + 0.4, y: "acf", fill: "#e67e22" })
  ]
})

Plot.plot({
  height: 350, marginLeft: 60,
  x: { label: "Lag", ticks: [5, 10, 15, 20, 25, 30, 35, 40, 45, 50] },
  y: { label: "Autocorrelation", grid: true, domain: [-0.15, 0.15] },
  marks: [
    Plot.ruleY([0], { stroke: "#ccc" }),
    Plot.ruleY([g22Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.ruleY([-g22Bartlett], { stroke: "#d62728", strokeDasharray: "4 3" }),
    Plot.rectY(g22ACFStdSq, { x1: d => d.lag - 0.4, x2: d => d.lag + 0.4, y: "acf", fill: "#e67e22" })
  ]
})

{
  const retStats = descStats(g22Sim.ret)
  const resid = g22Sim.ret.map((r, i) => r / Math.sqrt(g22Sim.gVar[i]))
  const rStats = descStats(resid)
  const lb15 = ljungBox(g22Data.map(d => d.stdSq), 15)

  return html`<table style="font-size:0.9rem;border-collapse:collapse;width:100%;">
    <thead>
      <tr style="border-bottom:2px solid #333;">
        <th style="text-align:left;padding:4px 8px;">Statistic</th>
        <th style="text-align:right;padding:4px 8px;">Returns</th>
        <th style="text-align:right;padding:4px 8px;">Std. residuals (z<sub>t</sub>)</th>
      </tr>
    </thead>
    <tbody>
      <tr><td style="padding:4px 8px;">Mean</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.mean, 6)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.mean, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Std. deviation</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.sd, 6)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.sd, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Skewness</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.skew, 4)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.skew, 4)}</td></tr>
      <tr><td style="padding:4px 8px;">Kurtosis</td><td style="text-align:right;padding:4px 8px;">${fmt(retStats.kurt, 2)}</td><td style="text-align:right;padding:4px 8px;">${fmt(rStats.kurt, 4)}</td></tr>
      <tr style="border-top:1px solid #ccc;">
        <td style="padding:4px 8px;">Ljung-Box(15) of R²/σ²</td><td style="text-align:right;padding:4px 8px;"></td><td style="text-align:right;padding:4px 8px;">${fmt(lb15, 2)}</td></tr>
    </tbody>
    <tfoot>
      <tr style="border-top:2px solid #333;">
        <td colspan="3" style="padding:4px 8px;font-size:0.82rem;color:#666;">
          Persistence = ${fmt(g22Pers, 2)}. Long-run vol = ${g22LRVol}%/day.
          ${!g22Stationary ? '<strong style="color:#d62728;">Model is non-stationary (persistence ≥ 1).</strong>' : ''}
          LB critical value (95%, K=15) = 25.0.
        </td>
      </tr>
    </tfoot>
  </table>`
}

4. Multi-period variance forecast

A key advantage of GARCH over EWMA is mean-reverting forecasts. The expected variance \(k\) days ahead under GARCH(1,1) is:

\[ E_t[\sigma^2_{t+k}] = V_L + (\alpha + \beta)^{k-1}\,(\sigma^2_{t+1} - V_L) \]

As \(k \to \infty\), the forecast converges to \(V_L\). Under EWMA (\(\alpha + \beta = 1\)), the forecast is flat: \(E_t[\sigma^2_{t+k}] = \sigma^2_{t+1}\) for all \(k\).

For multi-day VaR, we need the cumulative variance of the \(K\)-day return \(R_{t+1:t+K}\):

\[ \sigma^2_{t+1:t+K} = K\,V_L + \frac{1-(\alpha+\beta)^K}{1-(\alpha+\beta)}\,(\sigma^2_{t+1} - V_L) \]

Tip

How to experiment

Start with current volatility above the long-run level and watch the GARCH forecast curve downward while EWMA stays flat. Then try setting current volatility below long-run to see upward reversion. Increase persistence toward 1 to see the two models converge.

viewof fcPers = Inputs.range([0.90, 0.999], {
  label: "α + β (persistence)",
  step: 0.001,
  value: 0.970
})

viewof fcCurVol = Inputs.range([0.5, 4.0], {
  label: "Current daily volatility σ₀ (%)",
  step: 0.1,
  value: 2.5
})

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

viewof fcK = Inputs.range([10, 500], {
  label: "Forecast horizon K (days)",
  step: 10,
  value: 250
})

fcVL = (fcLRVol / 100) ** 2
fcV0 = (fcCurVol / 100) ** 2

// Compute forecast paths
fcData = {
  const data = []
  const p = fcPers, vl = fcVL, v0 = fcV0
  let cumGARCH = 0, cumEWMA = 0
  for (let k = 1; k <= fcK; k++) {
    const garchVar = vl + Math.pow(p, k - 1) * (v0 - vl)
    const ewmaVar = v0
    cumGARCH += garchVar
    cumEWMA += ewmaVar
    data.push({
      k,
      garchVol: Math.sqrt(garchVar),
      ewmaVol: Math.sqrt(ewmaVar),
      garchCum: cumGARCH,
      ewmaCum: cumEWMA
    })
  }
  return data
}

// Volatility term structure (annualized, Hull eq 8.15)
fcTermData = {
  const a = -Math.log(fcPers)
  const data = []
  for (let T = 1; T <= fcK; T++) {
    const avgVar = fcVL + (1 - Math.exp(-a * T)) / (a * T) * (fcV0 - fcVL)
    const annVol = Math.sqrt(252 * avgVar)
    data.push({ T, annVol })
  }
  return data
}

Plot.plot({
  height: 380, marginLeft: 65,
  x: { label: "Days ahead (k)", grid: false },
  y: { label: "Expected daily volatility", grid: true, tickFormat: pctFmt },
  marks: [
    Plot.ruleY([fcLRVol / 100], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
    Plot.line(fcData, { x: "k", y: "ewmaVol", stroke: "#d62728", strokeWidth: 2 }),
    Plot.line(fcData, { x: "k", y: "garchVol", stroke: "#2f71d5", strokeWidth: 2 })
  ]
})

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#2f71d5" stroke-width="2"/></svg> GARCH forecast</span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#d62728" stroke-width="2"/></svg> EWMA forecast (flat)</span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#999" stroke-width="1.5" stroke-dasharray="6 3"/></svg> Long-run volatility</span>
</div>
<p style="color:#666;font-size:0.85rem;">The GARCH forecast ${fcV0 > fcVL ? "declines" : "increases"} toward the long-run level. With persistence = ${fmt(fcPers, 3)}, the half-life of a shock is ${fmt(Math.log(2) / (-Math.log(fcPers)), 1)} days.</p>`

Plot.plot({
  height: 380, marginLeft: 65,
  x: { label: "Horizon K (days)", grid: false },
  y: { label: "Cumulative variance σ²(t+1:t+K)", grid: true },
  marks: [
    Plot.line(fcData, { x: "k", y: "ewmaCum", stroke: "#d62728", strokeWidth: 2 }),
    Plot.line(fcData, { x: "k", y: "garchCum", stroke: "#2f71d5", strokeWidth: 2 })
  ]
})

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#2f71d5" stroke-width="2"/></svg> GARCH cumulative</span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#d62728" stroke-width="2"/></svg> EWMA cumulative (= K × σ²₀)</span>
</div>
<p style="color:#666;font-size:0.85rem;">The cumulative variance determines multi-day VaR. When current volatility exceeds the long-run level, EWMA overstates the K-day risk relative to GARCH because it ignores mean reversion.</p>`

Plot.plot({
  height: 380, marginLeft: 65,
  x: { label: "Option maturity T (days)", grid: false },
  y: { label: "Annualized volatility σ(T)", grid: true, tickFormat: d => (d * 100).toFixed(0) + "%" },
  marks: [
    Plot.ruleY([Math.sqrt(252 * fcVL)], { stroke: "#999", strokeDasharray: "6 3", strokeWidth: 1.5 }),
    Plot.line(fcTermData, { x: "T", y: "annVol", stroke: "#2f71d5", strokeWidth: 2 })
  ]
})

html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#2f71d5" stroke-width="2"/></svg> GARCH term structure</span>
  <span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#999" stroke-width="1.5" stroke-dasharray="6 3"/></svg> Long-run annualized vol (${fmt(Math.sqrt(252 * fcVL) * 100, 1)}%)</span>
</div>
<p style="color:#666;font-size:0.85rem;">The volatility term structure for option pricing (Hull's equation 8.15), with a = −ln(α+β) = ${fmt(-Math.log(fcPers), 4)}. Current volatility ${fcV0 > fcVL ? "above" : "below"} long-run implies a ${fcV0 > fcVL ? "downward" : "upward"}-sloping term structure.</p>`

5. News impact function

The news impact function (NIF) shows how today’s standardized shock \(z_t\) affects tomorrow’s variance. Standard GARCH assumes a symmetric parabola, but the leverage effect (negative returns increase volatility more than positive returns) motivates asymmetric extensions (see Christoffersen 2012, sec. 5.1).

Model	Formula	NIF
GARCH	\(\sigma^2_{t+1} = \omega + \alpha\sigma^2_t z^2_t + \beta\sigma^2_t\)	\(z^2_t\) (symmetric)
NGARCH	\(\sigma^2_{t+1} = \omega + \alpha\sigma^2_t(z_t - \theta)^2 + \beta\sigma^2_t\)	\((z_t - \theta)^2\) (shifted)
GJR-GARCH	\(\sigma^2_{t+1} = \omega + \alpha R^2_t + \alpha\theta I_t R^2_t + \beta\sigma^2_t\)	\(z^2_t + \theta I_t z^2_t\) (kinked)

Tip

How to experiment

Select different models and adjust the leverage parameter \(\theta\). With NGARCH, a positive \(\theta\) shifts the minimum of the parabola to the right, so negative shocks (left of the minimum) produce larger variance increases. With GJR-GARCH, \(\theta > 0\) adds an extra kick for negative returns, creating a visible kink at \(z_t = 0\).

viewof nifModel = Inputs.radio(["GARCH", "NGARCH", "GJR-GARCH", "All models"], {
  label: "Model",
  value: "All models"
})

viewof nifAlpha = Inputs.range([0.01, 0.20], { label: "α", step: 0.01, value: 0.10 })

viewof nifBeta = Inputs.range([0.50, 0.95], { label: "β", step: 0.01, value: 0.85 })

viewof nifTheta = Inputs.range([0.0, 3.0], { label: "θ (leverage parameter)", step: 0.1, value: 1.0 })

viewof nifSigma = Inputs.range([0.5, 3.0], { label: "Current daily volatility σ_t (%)", step: 0.1, value: 1.5 })

nifSigma2 = (nifSigma / 100) ** 2

// Compute NIF data
nifData = {
  const a = nifAlpha, b = nifBeta, th = nifTheta, s2 = nifSigma2
  const pts = []

  // Long-run variances for omega calculation
  const omegaG = s2 * (1 - a - b)
  const omegaN = s2 * (1 - a * (1 + th * th) - b)
  const omegaGJR = s2 * (1 - a - 0.5 * a * th - b)  // approximate: E[I_t] = 0.5

  for (let z = -4; z <= 4; z += 0.05) {
    const R = z * Math.sqrt(s2)
    const I = z < 0 ? 1 : 0

    const garch = omegaG + a * R * R + b * s2
    const ngarch = omegaN + a * s2 * (z - th) ** 2 + b * s2
    const gjr = omegaGJR + a * R * R + a * th * I * R * R + b * s2

    pts.push({
      z,
      GARCH: Math.sqrt(garch) * 100,
      NGARCH: Math.sqrt(Math.max(0, ngarch)) * 100,
      "GJR-GARCH": Math.sqrt(Math.max(0, gjr)) * 100
    })
  }
  return pts
}

{
  const show = nifModel
  const models = show === "All models"
    ? ["GARCH", "NGARCH", "GJR-GARCH"]
    : [show]
  const colors = { "GARCH": "#888", "NGARCH": "#2f71d5", "GJR-GARCH": "#d62728" }

  const marks = []
  // Always show GARCH as reference if showing a specific asymmetric model
  if (show !== "All models" && show !== "GARCH") {
    marks.push(Plot.line(nifData, {
      x: "z", y: "GARCH", stroke: "#ccc", strokeWidth: 1.5, strokeDasharray: "6 3"
    }))
  }
  for (const m of models) {
    marks.push(Plot.line(nifData, {
      x: "z", y: m, stroke: colors[m], strokeWidth: 2
    }))
  }
  // Add vertical line at z=0
  marks.unshift(Plot.ruleX([0], { stroke: "#ddd" }))

  return Plot.plot({
    height: 400, marginLeft: 65,
    x: { label: "Standardized shock z_t", grid: false },
    y: { label: "Tomorrow's volatility σ(t+1) (%)", grid: true },
    marks
  })
}

{
  const show = nifModel
  const items = []
  if (show === "All models" || show === "GARCH")
    items.push(html`<span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#888" stroke-width="2"/></svg> GARCH (symmetric)</span>`)
  if (show !== "All models" && show !== "GARCH")
    items.push(html`<span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#ccc" stroke-width="1.5" stroke-dasharray="6 3"/></svg> GARCH (reference)</span>`)
  if (show === "All models" || show === "NGARCH")
    items.push(html`<span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#2f71d5" stroke-width="2"/></svg> NGARCH</span>`)
  if (show === "All models" || show === "GJR-GARCH")
    items.push(html`<span><svg width="24" height="10"><line x1="0" y1="5" x2="24" y2="5" stroke="#d62728" stroke-width="2"/></svg> GJR-GARCH</span>`)

  return html`<div style="display:flex;gap:18px;font-size:0.85rem;margin-top:-6px;flex-wrap:wrap;">${items}</div>`
}

html`<p style="color:#666;font-size:0.85rem;">The news impact function shows σ<sub>t+1</sub> as a function of the standardized shock z<sub>t</sub>, holding current volatility at ${nifSigma}%. GARCH is symmetric around z<sub>t</sub> = 0. NGARCH shifts the minimum to z<sub>t</sub> = θ = ${fmt(nifTheta, 1)}. GJR-GARCH has a kink at z<sub>t</sub> = 0, with steeper slope for negative shocks.</p>`

{
  const a = nifAlpha, b = nifBeta, th = nifTheta, s2 = nifSigma2
  const omegaG = s2 * (1 - a - b)
  const omegaN = s2 * (1 - a * (1 + th * th) - b)
  const omegaGJR = s2 * (1 - a - 0.5 * a * th - b)
  const shocks = [-3, -2, -1, 1, 2, 3]

  const rows = shocks.map(z => {
    const R = z * Math.sqrt(s2)
    const I = z < 0 ? 1 : 0
    const gV = Math.sqrt(omegaG + a * R * R + b * s2) * 100
    const nV = Math.sqrt(Math.max(0, omegaN + a * s2 * (z - th) ** 2 + b * s2)) * 100
    const gjrV = Math.sqrt(Math.max(0, omegaGJR + a * R * R + a * th * I * R * R + b * s2)) * 100
    return html`<tr>
      <td style="text-align:center;padding:4px 8px;${z < 0 ? 'color:#d62728;font-weight:500;' : ''}">${z > 0 ? '+' : ''}${z}σ</td>
      <td style="text-align:right;padding:4px 8px;">${fmt(gV, 3)}%</td>
      <td style="text-align:right;padding:4px 8px;">${fmt(nV, 3)}%</td>
      <td style="text-align:right;padding:4px 8px;">${fmt(gjrV, 3)}%</td>
    </tr>`
  })

  return html`<table style="font-size:0.9rem;border-collapse:collapse;width:100%;">
    <thead>
      <tr style="border-bottom:2px solid #333;">
        <th style="text-align:center;padding:4px 8px;">Shock z<sub>t</sub></th>
        <th style="text-align:right;padding:4px 8px;">GARCH σ<sub>t+1</sub></th>
        <th style="text-align:right;padding:4px 8px;">NGARCH σ<sub>t+1</sub></th>
        <th style="text-align:right;padding:4px 8px;">GJR σ<sub>t+1</sub></th>
      </tr>
    </thead>
    <tbody>${rows}</tbody>
    <tfoot>
      <tr style="border-top:2px solid #333;">
        <td colspan="4" style="padding:4px 8px;font-size:0.82rem;color:#666;">
          Current σ<sub>t</sub> = ${nifSigma}%. α = ${nifAlpha}, β = ${nifBeta}, θ = ${fmt(nifTheta, 1)}.
          In GARCH, the response to ±kσ is identical. In NGARCH and GJR, negative shocks produce higher volatility than positive shocks of equal magnitude.
        </td>
      </tr>
    </tfoot>
  </table>`
}

References

Bollerslev, Tim. 1986. “Generalized Autoregressive Conditional Heteroskedasticity.” Journal of Econometrics 31 (3): 307–27.

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

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