Multivariate FHS

Interactive exploration of Multivariate Filtered Historical Simulation with constant correlations for portfolio risk measurement

Multivariate Filtered Historical Simulation extends the univariate FHS approach to portfolios of multiple assets (Christoffersen 2012, chap. 8). The key idea: estimate individual GARCH models for each asset, extract standardized residuals, and then draw entire shock vectors from the same historical day to preserve the cross-asset correlation structure.

Under constant correlations, the procedure is:

1. Estimate GARCH for each asset \(i\) to obtain \(\hat{\sigma}_{i,t}\)
2. Compute standardized residuals \(\hat{z}_{i,t} = R_{i,t}/\hat{\sigma}_{i,t}\)
3. Store shock vectors \(\{\hat{z}_{t+1-\tau}\}_{\tau=1}^{m}\) from the same day, where \(m\) is the number of historical observations
4. Draw entire shock vectors with replacement
5. Simulate returns: \(\hat{R}_{j,t+k} = \hat{\sigma}_{j,t+k} \cdot \hat{z}_{j,\tau(k)}\), updating GARCH variances
6. Compute portfolio returns: \(\hat{R}^{PF}_{t+k} = \sum_j w_j \hat{R}_{j,t+k}\)
7. Calculate VaR and ES from the simulated distribution

Note

Why draw entire vectors? Drawing each asset’s shock independently would destroy the historical correlation structure, leading to systematic underestimation of portfolio tail risk. By drawing all assets’ shocks from the same historical day, we preserve whatever dependence (linear and nonlinear) existed in the data.

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, 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 stdTSample(rng, d) {
    const normal = boxMuller(rng)
    const chi2 = 2 * gammaSample(rng, d / 2)
    return normal / Math.sqrt(chi2 / d) * Math.sqrt((d - 2) / d)
  }
  return { mulberry32, boxMuller, stdTSample }
}

fmt = (x, d) => x === undefined || isNaN(x) ? "N/A" : x.toFixed(d)
pctFmt = x => (x * 100).toFixed(2) + "%"

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)
  }
}

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

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 === "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 === "area") swatchHTML = `<svg width="14" height="14"><rect width="14" height="14" fill="${d.color}" opacity="0.3"/></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
}

Same-day draws vs. independent draws

The critical implementation choice in multivariate FHS is whether to draw shocks for all assets from the same historical day or to draw them independently. This illustration shows the dramatic impact on portfolio risk measurement.

Tip

How to experiment

Start with a positive correlation (\(\rho\) around 0.5) and compare the distributions in the “Same-day vs. independent” tab. The independent draw distribution will have thinner tails because it destroys the correlation. Increase the correlation to amplify the difference. Also try increasing the horizon to see how the effect compounds over multiple days.

viewof mfRho = Inputs.range([-0.50, 0.90], {
  label: "ρ₁₂ (true correlation)",
  step: 0.01,
  value: 0.50
})

viewof mfW1 = Inputs.range([0.0, 1.0], {
  label: "w₁ (weight of asset 1)",
  step: 0.01,
  value: 0.60
})

viewof mfK = Inputs.range([1, 20], {
  label: "Horizon K (days)",
  step: 1,
  value: 1
})

viewof mfFH = Inputs.range([1000, 10000], {
  label: "FHS simulation paths",
  step: 1000,
  value: 5000
})

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

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

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

mutable mfSeed = 0

viewof mfReseed = {
  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 mfSeed++ }
  return btn
}

mfAlpha1 = 0.08
mfBeta1 = 0.88
mfVL1 = (0.015)**2
mfOmega1 = mfVL1 * (1 - mfAlpha1 - mfBeta1)

mfAlpha2 = 0.10
mfBeta2 = 0.85
mfVL2 = (0.020)**2
mfOmega2 = mfVL2 * (1 - mfAlpha2 - mfBeta2)

// Generate historical data and extract standardized residuals
mfHist = {
  mfSeed
  const rng = rng_utils.mulberry32(77 + mfSeed)
  const N = 1500
  const rho = mfRho
  const sq = Math.sqrt(Math.max(0, 1 - rho * rho))

  const R1 = new Array(N), R2 = new Array(N)
  const sig2_1 = new Array(N), sig2_2 = new Array(N)
  const z1 = new Array(N), z2 = new Array(N)

  sig2_1[0] = mfVL1; sig2_2[0] = mfVL2

  for (let t = 0; t < N; t++) {
    let zu1, zu2
    if (mfDist === "Normal") {
      zu1 = rng_utils.boxMuller(rng); zu2 = rng_utils.boxMuller(rng)
    } else {
      zu1 = rng_utils.stdTSample(rng, mfDf); zu2 = rng_utils.stdTSample(rng, mfDf)
    }
    const s1 = zu1, s2 = rho * zu1 + sq * zu2

    R1[t] = Math.sqrt(sig2_1[t]) * s1
    R2[t] = Math.sqrt(sig2_2[t]) * s2
    z1[t] = R1[t] / Math.sqrt(sig2_1[t])
    z2[t] = R2[t] / Math.sqrt(sig2_2[t])

    if (t < N - 1) {
      sig2_1[t + 1] = mfOmega1 + mfAlpha1 * R1[t]**2 + mfBeta1 * sig2_1[t]
      sig2_2[t + 1] = mfOmega2 + mfAlpha2 * R2[t]**2 + mfBeta2 * sig2_2[t]
    }
  }

  // Current volatilities (end of sample)
  const curSig1 = Math.sqrt(sig2_1[N - 1])
  const curSig2 = Math.sqrt(sig2_2[N - 1])

  return { z1, z2, R1, R2, sig2_1, sig2_2, curSig1, curSig2, N }
}

// Run multivariate FHS: same-day draws vs independent draws
mfResults = {
  const rng = rng_utils.mulberry32(200 + mfSeed)
  const { z1, z2, curSig1, curSig2, N } = mfHist
  const FH = mfFH, K = mfK
  const w1 = mfW1, w2 = 1 - mfW1
  const p = mfP / 100

  const pfRetSame = new Array(FH)
  const pfRetIndep = new Array(FH)

  for (let i = 0; i < FH; i++) {
    let cumSame = 0, cumIndep = 0
    let s2_1 = curSig1 * curSig1
    let s2_2 = curSig2 * curSig2
    let s2_1i = s2_1, s2_2i = s2_2

    for (let k = 0; k < K; k++) {
      // Same-day draw
      const idx = Math.floor(rng() * N)
      const zs1 = z1[idx], zs2 = z2[idx]
      const rs1 = Math.sqrt(s2_1) * zs1
      const rs2 = Math.sqrt(s2_2) * zs2
      cumSame += w1 * rs1 + w2 * rs2
      s2_1 = mfOmega1 + mfAlpha1 * rs1**2 + mfBeta1 * s2_1
      s2_2 = mfOmega2 + mfAlpha2 * rs2**2 + mfBeta2 * s2_2

      // Independent draws
      const idx1 = Math.floor(rng() * N)
      const idx2 = Math.floor(rng() * N)
      const zi1 = z1[idx1], zi2 = z2[idx2]
      const ri1 = Math.sqrt(s2_1i) * zi1
      const ri2 = Math.sqrt(s2_2i) * zi2
      cumIndep += w1 * ri1 + w2 * ri2
      s2_1i = mfOmega1 + mfAlpha1 * ri1**2 + mfBeta1 * s2_1i
      s2_2i = mfOmega2 + mfAlpha2 * ri2**2 + mfBeta2 * s2_2i
    }

    pfRetSame[i] = cumSame
    pfRetIndep[i] = cumIndep
  }

  // Sort for percentile computation
  const sortedSame = pfRetSame.slice().sort((a, b) => a - b)
  const sortedIndep = pfRetIndep.slice().sort((a, b) => a - b)

  const posIdx = Math.max(0, Math.ceil(FH * p) - 1)
  const varSame = -sortedSame[posIdx]
  const varIndep = -sortedIndep[posIdx]

  // ES
  const tailSame = sortedSame.slice(0, posIdx + 1)
  const tailIndep = sortedIndep.slice(0, posIdx + 1)
  const esSame = -tailSame.reduce((a, b) => a + b, 0) / tailSame.length
  const esIndep = -tailIndep.reduce((a, b) => a + b, 0) / tailIndep.length

  // Normal VaR for comparison
  const pfSig = Math.sqrt(
    w1*w1*curSig1*curSig1 + w2*w2*curSig2*curSig2 + 2*w1*w2*curSig1*curSig2*mfRho
  )
  const normVaR = pfSig * Math.abs(qnorm(p)) * Math.sqrt(K)
  const normES = pfSig * normalPDF(qnorm(p)) / p * Math.sqrt(K)

  return {
    pfRetSame, pfRetIndep, sortedSame, sortedIndep,
    varSame, varIndep, esSame, esIndep,
    normVaR, normES, pfSig, p, K, FH
  }
}

Same-day vs. independent draws

viewof mfHistLegend = legend([
  { key: "same", label: "Same-day draws (correct)", color: "#2f71d5", type: "area" },
  { key: "indep", label: "Independent draws (incorrect)", color: "#e67e22", type: "area" }
])

{
  const h = mfHistLegend
  const marks = [Plot.ruleY([0])]

  if (!h.has("same"))
    marks.push(Plot.rectY(mfResults.pfRetSame, Plot.binX({ y: "count" }, {
      x: d => d, fill: "#2f71d5", fillOpacity: 0.4, thresholds: 80
    })))
  if (!h.has("indep"))
    marks.push(Plot.rectY(mfResults.pfRetIndep, Plot.binX({ y: "count" }, {
      x: d => d, fill: "#e67e22", fillOpacity: 0.3, thresholds: 80
    })))

  // VaR lines
  if (!h.has("same"))
    marks.push(Plot.ruleX([-mfResults.varSame], { stroke: "#2f71d5", strokeWidth: 2 }))
  if (!h.has("indep"))
    marks.push(Plot.ruleX([-mfResults.varIndep], { stroke: "#e67e22", strokeWidth: 2, strokeDasharray: "6 3" }))

  return Plot.plot({
    height: 400, marginLeft: 60,
    x: { label: `${mfK}-day portfolio return`, grid: false, tickFormat: pctFmt },
    y: { label: "Count", grid: true },
    marks
  })
}

{
  const diff = ((mfResults.varSame - mfResults.varIndep) / mfResults.varSame * 100)
  return html`<p style="color:#666;font-size:0.85rem;">Distribution of simulated ${mfK}-day portfolio returns. The <span style="color:#2f71d5;font-weight:700;">blue distribution</span> draws entire shock vectors from the same day (correct), preserving the correlation ρ = ${fmt(mfRho, 2)}. The <span style="color:#e67e22;font-weight:700;">orange distribution</span> draws each asset independently (incorrect), destroying the correlation. Vertical lines mark the respective VaR levels. The independent approach ${mfRho > 0 ? "underestimates" : "overestimates"} risk by ${fmt(Math.abs(diff), 1)}%.</p>`
}

VaR term structure

// Compute VaR across horizons
mfTermStructure = {
  const rng = rng_utils.mulberry32(300 + mfSeed)
  const { z1, z2, curSig1, curSig2, N } = mfHist
  const FH = Math.min(mfFH, 3000)
  const w1 = mfW1, w2 = 1 - mfW1
  const p = mfP / 100
  const maxK = 20

  const result = []
  for (let K = 1; K <= maxK; K++) {
    const pfRet = new Array(FH)
    for (let i = 0; i < FH; i++) {
      let cum = 0, s2_1 = curSig1**2, s2_2 = curSig2**2
      for (let k = 0; k < K; k++) {
        const idx = Math.floor(rng() * N)
        const rs1 = Math.sqrt(s2_1) * z1[idx]
        const rs2 = Math.sqrt(s2_2) * z2[idx]
        cum += w1 * rs1 + w2 * rs2
        s2_1 = mfOmega1 + mfAlpha1 * rs1**2 + mfBeta1 * s2_1
        s2_2 = mfOmega2 + mfAlpha2 * rs2**2 + mfBeta2 * s2_2
      }
      pfRet[i] = cum
    }
    pfRet.sort((a, b) => a - b)
    const posIdx = Math.max(0, Math.ceil(FH * p) - 1)
    const fhsVaR = -pfRet[posIdx]
    const tail = pfRet.slice(0, posIdx + 1)
    const fhsES = -tail.reduce((a, b) => a + b, 0) / tail.length

    // Normal scaling
    const pfSig = Math.sqrt(
      w1*w1*curSig1**2 + w2*w2*curSig2**2 + 2*w1*w2*curSig1*curSig2*mfRho
    )
    const normVaR = pfSig * Math.abs(qnorm(p)) * Math.sqrt(K)
    const normES = pfSig * normalPDF(qnorm(p)) / p * Math.sqrt(K)

    result.push({
      K,
      fhsVaR: fhsVaR * 100, fhsES: fhsES * 100,
      normVaR: normVaR * 100, normES: normES * 100,
      fhsVaRPerDay: fhsVaR / Math.sqrt(K) * 100,
      normVaRPerDay: normVaR / Math.sqrt(K) * 100
    })
  }
  return result
}

viewof mfTSLegend = legend([
  { key: "fhsvar", label: "FHS VaR", color: "#2f71d5", type: "line" },
  { key: "fhses", label: "FHS ES", color: "#d62728", type: "line" },
  { key: "normvar", label: "Normal VaR (√K scaling)", color: "#2f71d5", type: "dashed" },
  { key: "normes", label: "Normal ES (√K scaling)", color: "#d62728", type: "dashed" }
])

{
  const h = mfTSLegend
  const marks = []
  if (!h.has("normvar"))
    marks.push(Plot.line(mfTermStructure, { x: "K", y: "normVaR", stroke: "#2f71d5", strokeWidth: 1.5, strokeDasharray: "6 3" }))
  if (!h.has("normes"))
    marks.push(Plot.line(mfTermStructure, { x: "K", y: "normES", stroke: "#d62728", strokeWidth: 1.5, strokeDasharray: "6 3" }))
  if (!h.has("fhsvar"))
    marks.push(Plot.line(mfTermStructure, { x: "K", y: "fhsVaR", stroke: "#2f71d5", strokeWidth: 2.5 }))
  if (!h.has("fhses"))
    marks.push(Plot.line(mfTermStructure, { x: "K", y: "fhsES", stroke: "#d62728", strokeWidth: 2.5 }))

  // Highlight current K
  marks.push(Plot.ruleX([mfK], { stroke: "#999", strokeDasharray: "2 3" }))

  return Plot.plot({
    height: 400, marginLeft: 60,
    x: { label: "Horizon K (days)", grid: false, ticks: [1, 5, 10, 15, 20] },
    y: { label: "Risk measure (% of portfolio)", grid: true },
    marks
  })
}

html`<p style="color:#666;font-size:0.85rem;">Term structure of portfolio VaR and ES from multivariate FHS (solid) vs. normal √K scaling (dashed). With fat-tailed shocks, FHS produces higher risk estimates, especially at shorter horizons. The vertical gray line marks the currently selected horizon K = ${mfK}.</p>`

Summary

{
  const r = mfResults
  return html`<table class="table" style="width:100%;">
<thead><tr>
  <th>Measure</th>
  <th style="color:#2f71d5;">Same-day draws</th>
  <th style="color:#e67e22;">Independent draws</th>
  <th>Normal</th>
</tr></thead>
<tbody>
<tr><td style="font-weight:500;">${mfK}-day VaR (${fmt(mfP,1)}%)</td>
    <td style="font-weight:700;">${fmt(r.varSame * 100, 3)}%</td>
    <td>${fmt(r.varIndep * 100, 3)}%</td>
    <td>${fmt(r.normVaR * 100, 3)}%</td></tr>
<tr><td style="font-weight:500;">${mfK}-day ES (${fmt(mfP,1)}%)</td>
    <td style="font-weight:700;">${fmt(r.esSame * 100, 3)}%</td>
    <td>${fmt(r.esIndep * 100, 3)}%</td>
    <td>${fmt(r.normES * 100, 3)}%</td></tr>
<tr><td style="font-weight:500;">VaR underestimation (independent)</td>
    <td colspan="3" style="color:#d62728;font-weight:700;">${fmt((r.varSame - r.varIndep) / r.varSame * 100, 1)}%</td></tr>
<tr><td style="font-weight:500;">ES underestimation (independent)</td>
    <td colspan="3" style="color:#d62728;font-weight:700;">${fmt((r.esSame - r.esIndep) / r.esSame * 100, 1)}%</td></tr>
</tbody></table>
<p style="color:#666;font-size:0.85rem;">Current portfolio: w₁ = ${fmt(mfW1, 2)}, w₂ = ${fmt(1-mfW1, 2)}. Current volatilities: σ₁ = ${fmt(mfHist.curSig1 * 100, 3)}%, σ₂ = ${fmt(mfHist.curSig2 * 100, 3)}%. Correlation: ρ = ${fmt(mfRho, 2)}. Simulated with ${mfFH} paths over ${mfK} days.</p>`
}

References

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