Four Option VaR Methods Compared

Delta, gamma (Cornish-Fisher and simulated), and full-valuation VaR for a single-option portfolio, illustrated with a common set of Monte Carlo draws

For a position of \(m\) contracts (negative for a short position) in a single option on one underlying, the dollar change in portfolio value over a \(K\)-trading-day horizon is (Christoffersen 2012, chap. 11)

\[ \Delta P = m\,\bigl(c(S_{t+\tau}, \tilde T - \tau) - c(S_t, \tilde T)\bigr), \]

where \(c(S, T)\) is the option price as a function of spot \(S\) and time to maturity \(T\), \(S_t\) is today’s spot, \(\tilde T\) is today’s time to maturity in calendar days, and \(\tau\) is the risk horizon in calendar days (with \(\tau = K\times 365/252\) when \(K\) is measured in trading days). Four common VaR estimators at tail probability \(p\) are:

Note

Convention on \(\tau\). This page uses the exact conversion \(\tau = K\times 365/252 \approx 1.4484\,K\), which preserves the annualized time. Christoffersen (2012) rounds this to \(\tau = 14\) calendar days for \(K = 10\) trading days (a ratio of \(1.4\)). For short horizons the difference is numerically small but non-zero, so VaR values here will not match the textbook to the last decimal.

1. Delta-based (analytical). Linearize \(\Delta P\) around \(S_t\) and assume a normal return: \[ \text{VaR}^{\,p} = |m\,\delta|\,S_t\,\sigma_{\text{tr}}\sqrt{K}\,\lvert\Phi_p^{-1}\rvert, \] where \(\delta = \partial c/\partial S\) is today’s option delta, \(\sigma_{\text{tr}}\) is the daily trading-day volatility of the underlying, and \(\Phi_p^{-1}\) is the \(p\)-quantile of the standard normal distribution.

2. Gamma-based, Cornish–Fisher. Add the gamma term \(\gamma = \partial^{2} c/\partial S^{2}\), compute the first three moments of the quadratic approximation analytically, and apply the Cornish–Fisher expansion: \[ \text{VaR}^{\,p} = -\mu_{\Delta P} - \Bigl(\Phi_p^{-1} + \tfrac{1}{6}\!\bigl((\Phi_p^{-1})^{2} - 1\bigr)\,\zeta_{1,\Delta P}\Bigr)\,\sigma_{\Delta P}, \] where \(\mu_{\Delta P}\), \(\sigma_{\Delta P}\), and \(\zeta_{1,\Delta P}\) are the mean, standard deviation, and skewness of \(\Delta P\) under the quadratic approximation.

3. Gamma-based, simulated. Draw \(R_h \sim \mathcal N(0, K\sigma_{\text{tr}}^{2})\) for \(h = 1, \dots, \text{MC}\) and compute \(\widehat{\Delta P}_h = m\,\delta\,S_t R_h + \tfrac{1}{2}\,m\,\gamma\,S_t^{2} R_h^{2}\); take the empirical \(p\)-quantile.

4. Full valuation. Use the same \(R_h\), set \(\hat S_h = S_t\,e^{R_h}\), and reprice the option at the reduced maturity: \(\widehat{\Delta P}_h = m\,\bigl(c(\hat S_h, \tilde T - \tau) - c(S_t, \tilde T)\bigr)\).

The first three use local greeks; only full valuation handles the true non-linearity, including time decay.

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

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

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

bsm = (S, X, T, r, q, sig, type) => {
  if (T <= 0) {
    const intrinsic = type === "call" ? Math.max(S - X, 0) : Math.max(X - S, 0)
    return { price: intrinsic, delta: type === "call" ? (S > X ? 1 : 0) : (S < X ? -1 : 0), gamma: 0 }
  }
  const sqrtT = Math.sqrt(T)
  const d1 = (Math.log(S / X) + (r - q + 0.5 * sig * sig) * T) / (sig * sqrtT)
  const d2 = d1 - sig * sqrtT
  const eqT = Math.exp(-q * T), erT = Math.exp(-r * T)
  const Nd1 = normalCDF(d1), Nd2 = normalCDF(d2), pdf = normalPDF(d1)
  let price, delta
  if (type === "call") {
    price = S * eqT * Nd1 - X * erT * Nd2
    delta = eqT * Nd1
  } else {
    price = X * erT * normalCDF(-d2) - S * eqT * normalCDF(-d1)
    delta = eqT * (Nd1 - 1)
  }
  const gamma = eqT * pdf / (S * sig * sqrtT)
  return { price, delta, gamma }
}

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

boxMullerPair = (u1, u2) => {
  const u = Math.max(u1, 1e-12)
  const mag = Math.sqrt(-2 * Math.log(u))
  return [mag * Math.cos(2 * Math.PI * u2), mag * Math.sin(2 * Math.PI * u2)]
}

fmt = (x, d) => (x === undefined || !isFinite(x) || isNaN(x)) ? "N/A" : x.toFixed(d)

fmtMoney = x => (!isFinite(x) || isNaN(x)) ? "N/A" : (x < 0 ? "-$" : "$") + Math.abs(x).toFixed(2)

Setting up the portfolio

Tip

How to experiment

The defaults are a short one-contract at-the-money call on an 80-unit underlying. The short gamma pulls the true P&L distribution sharply negative, so the delta-based VaR understates risk. Move the option deep in-the-money and all four VaRs collapse onto each other: linearity returns. Flip to a long call and the linear approximation instead overstates downside risk, because gamma is now on your side. Change the horizon or volatility to see how the gap between methods scales.

viewof vmSt = Inputs.range([20, 200], { label: "Sₜ current spot", step: 0.5, value: 80 })

viewof vmX = Inputs.range([20, 200], { label: "X strike", step: 0.5, value: 80 })

viewof vmSig = Inputs.range([5, 80], { label: "σ annual volatility (%)", step: 0.5, value: 28 })

viewof vmDays = Inputs.range([10, 365], { label: "T̃ maturity (calendar days)", step: 1, value: 45 })

viewof vmR = Inputs.range([0, 8], { label: "r risk-free rate (%/year)", step: 0.1, value: 2.0 })

viewof vmQ = Inputs.range([0, 8], { label: "q dividend yield (%/year)", step: 0.1, value: 1.0 })

viewof vmType = Inputs.radio(["call", "put"], { label: "Option type", value: "call" })

viewof vmPosition = Inputs.range([-5, 5], { label: "m position (contracts; negative = short)", step: 0.5, value: -1 })

viewof vmK = Inputs.range([1, 30], { label: "K horizon (trading days)", step: 1, value: 10 })

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

viewof vmMC = Inputs.range([500, 20000], { label: "Monte Carlo draws MC", step: 500, value: 5000 })

mutable vmSeed = 12345

viewof vmReseed = {
  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 random draws</button>`
  btn.onclick = () => { mutable vmSeed++ }
  return btn
}

vmGreeks = {
  const S = vmSt, X = vmX
  const sig = vmSig / 100
  const T = vmDays / 365
  const r = vmR / 100, q = vmQ / 100
  const b = bsm(S, X, T, r, q, sig, vmType)
  return { S, X, T, sig, r, q, type: vmType, price: b.price, delta: b.delta, gamma: b.gamma }
}

vmSim = {
  const g = vmGreeks
  const m = vmPosition
  const K = vmK
  const MC = vmMC
  const p = vmTail / 100

  // Trading-day vol: σ_ann uses 252 trading days
  const sigTr = g.sig / Math.sqrt(252)
  const sigK = sigTr * Math.sqrt(K)

  // Generate N(0,1) draws via Box-Muller on mulberry32
  const rng = mulberry32(vmSeed)
  const z = []
  for (let i = 0; i < MC; i += 2) {
    const [a, b] = boxMullerPair(rng(), rng())
    z.push(a)
    if (i + 1 < MC) z.push(b)
  }

  // K-day return draws (arithmetic, used both as R for the quadratic approx and as the log return for full valuation)
  const R = z.map(v => v * sigK)

  // Reduced maturity for full valuation: τ = K × (365/252) calendar days
  const tauCal = K * 365 / 252
  const Tnew = Math.max((g.T * 365 - tauCal) / 365, 0.001 / 365)

  // Per-draw ΔP under the three simulation-based methods
  const dPDelta = new Array(MC)
  const dPGamma = new Array(MC)
  const dPFull  = new Array(MC)
  const cMkt = g.price
  const delta_e = m * g.delta
  const gamma_e = m * g.gamma

  for (let i = 0; i < MC; i++) {
    const Ri = R[i]
    dPDelta[i] = delta_e * g.S * Ri
    dPGamma[i] = delta_e * g.S * Ri + 0.5 * gamma_e * g.S * g.S * Ri * Ri
    const Snew = g.S * Math.exp(Ri)
    const cNew = bsm(Snew, g.X, Tnew, g.r, g.q, g.sig, g.type).price
    dPFull[i] = m * (cNew - cMkt)
  }

  // Analytical delta VaR
  const zp = Math.abs(qnorm(p))
  const varDeltaAnalytic = Math.abs(m * g.delta) * g.S * sigTr * Math.sqrt(K) * zp

  // Cornish-Fisher moments for gamma VaR
  const sigmaK2 = K * sigTr * sigTr
  const mu  = 0.5 * gamma_e * g.S * g.S * sigmaK2
  const var_ = delta_e * delta_e * g.S * g.S * sigmaK2
             + 0.5 * gamma_e * gamma_e * Math.pow(g.S, 4) * sigmaK2 * sigmaK2
  const sigma = Math.sqrt(Math.max(var_, 1e-30))
  const term1 = (9/2) * delta_e * delta_e * gamma_e * Math.pow(g.S, 4) * sigmaK2 * sigmaK2
  const term2 = (15/8) * Math.pow(gamma_e, 3) * Math.pow(g.S, 6) * Math.pow(sigmaK2, 3)
  const term3 = -3 * (delta_e * delta_e * g.S * g.S * sigmaK2
                    + 0.75 * gamma_e * gamma_e * Math.pow(g.S, 4) * sigmaK2 * sigmaK2) * mu
  const term4 = 2 * Math.pow(mu, 3)
  const skew = (term1 + term2 + term3 + term4) / Math.pow(sigma, 3)

  const zpSigned = qnorm(p)
  const varGammaCF = -mu - (zpSigned + (1/6) * (zpSigned * zpSigned - 1) * skew) * sigma

  // Empirical percentile for the three simulated distributions
  const percentile = (arr, p) => {
    const sorted = arr.slice().sort((x, y) => x - y)
    const idx = p * (sorted.length - 1)
    const lo = Math.floor(idx), hi = Math.ceil(idx)
    if (lo === hi) return sorted[lo]
    return sorted[lo] * (hi - idx) + sorted[hi] * (idx - lo)
  }
  const varDeltaSim = -percentile(dPDelta, p)
  const varGammaSim = -percentile(dPGamma, p)
  const varFullSim  = -percentile(dPFull,  p)

  // ES for full valuation (empirical)
  const cutoff = percentile(dPFull, p)
  const tailLosses = dPFull.filter(x => x <= cutoff)
  const esFullSim = tailLosses.length > 0
    ? -tailLosses.reduce((a, b) => a + b, 0) / tailLosses.length
    : NaN

  return {
    m, K, p, MC, sigTr, sigK, Tnew,
    R, dPDelta, dPGamma, dPFull,
    varDeltaAnalytic, varGammaCF,
    varDeltaSim, varGammaSim, varFullSim,
    esFullSim,
    mu, sigma, skew
  }
}

// Shared bins across methods for fair histogram comparison
vmBins = {
  const all = [...vmSim.dPDelta, ...vmSim.dPGamma, ...vmSim.dPFull]
  const lo = Math.min(...all)
  const hi = Math.max(...all)
  const n = 60
  const w = (hi - lo) / n
  const bin = (arr) => {
    const b = new Array(n).fill(0).map((_, i) => ({
      xMid: lo + (i + 0.5) * w, x0: lo + i * w, x1: lo + (i + 1) * w, count: 0
    }))
    for (const v of arr) {
      let i = Math.floor((v - lo) / w)
      if (i === n) i = n - 1
      if (i >= 0 && i < n) b[i].count++
    }
    return b
  }
  return {
    lo, hi,
    delta: bin(vmSim.dPDelta).map(d => ({ ...d, method: "Delta (sim)" })),
    gamma: bin(vmSim.dPGamma).map(d => ({ ...d, method: "Gamma (sim)" })),
    full:  bin(vmSim.dPFull ).map(d => ({ ...d, method: "Full valuation" }))
  }
}

P&L distributions

{
  const all = [...vmBins.delta, ...vmBins.gamma, ...vmBins.full]
  const methods = ["Delta (sim)", "Gamma (sim)", "Full valuation"]
  const colors = { "Delta (sim)": "#e67e22", "Gamma (sim)": "#2e7d32", "Full valuation": "#2f71d5" }

  const vaRLines = [
    { method: "Delta (sim)",    x: -vmSim.varDeltaSim, v: vmSim.varDeltaSim },
    { method: "Gamma (sim)",    x: -vmSim.varGammaSim, v: vmSim.varGammaSim },
    { method: "Full valuation", x: -vmSim.varFullSim,  v: vmSim.varFullSim  }
  ]

  return Plot.plot({
    height: 520, marginLeft: 60, marginRight: 130, marginBottom: 35,
    x: { label: "Portfolio ΔP ($)", grid: false },
    y: { label: "Frequency", grid: true },
    color: { domain: methods, range: methods.map(m => colors[m]) },
    fy: { domain: methods, label: null },
    marks: [
      Plot.rectY(all, {
        x1: "x0", x2: "x1", y: "count", fy: "method",
        fill: d => colors[d.method], fillOpacity: 0.75
      }),
      Plot.ruleY([0]),
      Plot.ruleX(vaRLines, {
        x: "x", fy: "method",
        stroke: "#c0392b", strokeWidth: 2, strokeDasharray: "4 3"
      }),
      Plot.text(vaRLines, {
        x: "x", fy: "method",
        dx: 6, dy: 14, fontSize: 10, fontWeight: 600, fill: "#c0392b",
        text: d => `VaR = ${fmtMoney(d.v)}`,
        textAnchor: "start", frameAnchor: "top-left"
      })
    ]
  })
}

html`<p style="color:#666;font-size:0.85rem;">Three P&L distributions built from the <em>same</em> ${fmt(vmSim.MC, 0)} return draws: the linear (delta) projection, the quadratic (gamma) projection, and the exact (full-valuation) reprice. Dashed red lines mark the empirical VaR at the ${fmt(vmSim.p*100, 1)}% tail. The delta distribution is symmetric by construction; gamma introduces skewness; full valuation is the benchmark.</p>`

VaR comparison

{
  const bars = [
    { method: "Delta (analytic)", value: vmSim.varDeltaAnalytic, idx: 0 },
    { method: "Delta (simulated)", value: vmSim.varDeltaSim, idx: 1 },
    { method: "Gamma (Cornish-Fisher)", value: vmSim.varGammaCF, idx: 2 },
    { method: "Gamma (simulated)", value: vmSim.varGammaSim, idx: 3 },
    { method: "Full valuation", value: vmSim.varFullSim, idx: 4 }
  ]
  const colors = ["#f39c12", "#e67e22", "#7cb342", "#2e7d32", "#2f71d5"]

  return Plot.plot({
    height: 360, marginLeft: 60, marginRight: 20, marginBottom: 80,
    x: { label: null, domain: bars.map(b => b.method), tickRotate: -20 },
    y: { label: `${vmSim.K}-day ${fmt(vmSim.p*100, 1)}% VaR ($)`, grid: true },
    marks: [
      Plot.ruleY([0]),
      Plot.barY(bars, { x: "method", y: "value", fill: (d, i) => colors[d.idx], fillOpacity: 0.85 }),
      Plot.text(bars, {
        x: "method", y: "value",
        text: d => fmtMoney(d.value),
        dy: -6, fontSize: 11, fontWeight: 600
      }),
      Plot.ruleY([vmSim.varFullSim], { stroke: "#2f71d5", strokeDasharray: "4 3", strokeOpacity: 0.5 })
    ]
  })
}

{
  const v = vmSim
  const full = v.varFullSim
  const ratio = x => full > 0 ? 100 * x / full : NaN
  return html`<p style="color:#666;font-size:0.85rem;">Full valuation (blue) is the benchmark. Delta-analytic / delta-sim converge as MC grows. Ratios to full valuation: delta-analytic ${fmt(ratio(v.varDeltaAnalytic), 0)}%, delta-sim ${fmt(ratio(v.varDeltaSim), 0)}%, gamma-CF ${fmt(ratio(v.varGammaCF), 0)}%, gamma-sim ${fmt(ratio(v.varGammaSim), 0)}%.</p>`
}

Moments and summary

{
  const g = vmGreeks, v = vmSim
  const el = html`<table class="table" style="width:100%;">
    <thead><tr><th>Quantity</th><th>Value</th></tr></thead>
    <tbody>
      <tr><td>Position \\(m\\)</td><td>${v.m}</td></tr>
      <tr><td>Option price at \\(S_t\\)</td><td>${fmt(g.price, 4)}</td></tr>
      <tr><td>Delta \\(\\delta\\) (per contract)</td><td>${fmt(g.delta, 4)}</td></tr>
      <tr><td>Gamma \\(\\gamma\\) (per contract)</td><td>${fmt(g.gamma, 6)}</td></tr>
      <tr><td>Portfolio \\(m\\,\\delta\\)</td><td>${fmt(v.m * g.delta, 4)}</td></tr>
      <tr><td>Portfolio \\(m\\,\\gamma\\)</td><td>${fmt(v.m * g.gamma, 6)}</td></tr>
      <tr><td colspan="2" style="background:#f5f5f5;font-weight:600;">Moments of \\(\\Delta P\\) (${v.K}-day, quadratic approximation)</td></tr>
      <tr><td>Mean \\(\\mu\\)</td><td>${fmtMoney(v.mu)}</td></tr>
      <tr><td>Std. dev. \\(\\sigma\\)</td><td>${fmtMoney(v.sigma)}</td></tr>
      <tr><td>Skewness \\(\\zeta_1\\)</td><td>${fmt(v.skew, 3)}</td></tr>
      <tr><td colspan="2" style="background:#f5f5f5;font-weight:600;">VaR estimates (${fmt(v.p*100, 1)}% tail, ${v.K}-day)</td></tr>
      <tr><td>Delta (analytical)</td><td>${fmtMoney(v.varDeltaAnalytic)}</td></tr>
      <tr><td>Delta (simulated)</td><td>${fmtMoney(v.varDeltaSim)}</td></tr>
      <tr><td>Gamma (Cornish-Fisher)</td><td>${fmtMoney(v.varGammaCF)}</td></tr>
      <tr><td>Gamma (simulated)</td><td>${fmtMoney(v.varGammaSim)}</td></tr>
      <tr><td>Full valuation</td><td><strong>${fmtMoney(v.varFullSim)}</strong></td></tr>
      <tr><td>Full valuation ES</td><td>${fmtMoney(v.esFullSim)}</td></tr>
    </tbody></table>`
  if (window.MathJax && MathJax.typesetPromise) MathJax.typesetPromise([el])
  return el
}

Note

Why can gamma-simulated and full-valuation disagree? Even though both use the same \(R_h\) draws, gamma-sim truncates the Taylor series at second order and ignores time decay. Full valuation reprices the option with reduced maturity \(\tilde T - \tau\), so it captures theta as well as higher-order curvature terms. The two methods converge as \(K \to 0\) and as the option moves away from the money.

References

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