Realized vs implied volatility; volatility risk premium; VIX and VSTOXX.

Two different volatilities

The word “volatility” refers to two things that should not be confused:

• Realized volatility is a statistical measure of the past: the annualized standard deviation of the returns of a period that has already happened. With daily log returns $r_t$ and $n$ sessions, $\sigma_{\text{real}} = \sqrt{(252/n)\sum r_t^2}$ (without subtracting the mean, as in variance contracts). There are better estimators using highs and lows (the site's “Volatility Estimators” research compares them), but the close-to-close one is the one derivatives settle on.
• Implied volatility is a price: what the options market charges today for the volatility of a future period (lessons 3.4.4 and 3.4.5). It is not a neutral forecast, because it includes what buyers pay to protect themselves.

The difference between today's implied and the realized that comes later is the volatility risk premium. If it is positive on average, selling volatility makes money on average; the question is in exchange for what risk. Measuring it requires a constant-tenor implied volatility, and the market standard is Cboe's VIX: the 30-day implied volatility of the S&P 500 computed with the whole option chain, not with a single option.

flowchart TB
  O[OTM options] --> V[Variance at T1 and T2]
  V --> I[Interpolate to 30 d]
  I --> X[VIX = 100 sqrt]
  X --> P{Compare}
  R[30-day realized] --> P
  P --> VRP[Risk premium]

From the chain to the index: Cboe's formula

The VIX was born in 1993 on ATM options on the S&P 100. In 2003 Cboe, with Goldman Sachs, rebuilt it on S&P 500 options with the current formula, which does not depend on any pricing model11Note 1 · History and current rules: Cboe, Volatility Index Methodology: Cboe Volatility Index (version 6.0) and Cboe Volatility Index Mathematics Methodology; monthly version with at least one week to go: 2009 white paper. Links in the lesson resources.. The idea comes from the theory of variance swaps (Demeterfi, Derman, Kamal and Zou, 1999; Carr and Madan, 1998): a portfolio of out-of-the-money options with weights inversely proportional to the square of the strike replicates realized variance, whatever the level of the index. Its price, therefore, is the expected variance under the pricing measure. For an expiry $T$:

Each term has its rule in Cboe's methodology:

Then the two expiries are interpolated to 30 days in total variance and converted to volatility:

The real VIX today uses options with more than 23 and fewer than 37 days, monthly and weekly, and counts time in minutes. The 2009 white paper, with only monthly expiries, took the first two expiries with at least one week to go. That monthly version is the one used by the VORX, the VIX-like index of the course data on the fictitious ORX index (vol3_indices.csv): the first two expiries with at least 8 days, tenor in calendar days and the rate linearly interpolated on the curve in opciones3_indice.csv. Its known answer is that it can be recomputed exactly from the published files.

One day's VORX, step by step

import math
import numpy as np
import pandas as pd

cad = pd.read_parquet("opciones3_cadena.parquet")
ind = pd.read_csv("opciones3_indice.csv", parse_dates=["fecha"]) \
    .set_index("fecha")
vx = pd.read_csv("vol3_indices.csv", parse_dates=["fecha"]) \
    .set_index("fecha")

ancha = cad.pivot_table(index=["fecha", "vencimiento", "strike"],
                        columns="tipo", values=["bid", "ask"])
ancha.columns = [f"{a}_{t}" for a, t in ancha.columns]
ancha = ancha.dropna().reset_index()   # strikes with call and put

PLAZOS = np.array([1 / 12, 0.25, 0.5, 1.0, 2.0])
CURVA = ["tipo_1m", "tipo_3m", "tipo_6m", "tipo_1a", "tipo_2a"]


def tipo(fecha, T):
    """Continuous rate at tenor T: linear interpolation of the curve."""
    return float(np.interp(T, PLAZOS, ind.loc[fecha, CURVA])) / 100


hoy = pd.Timestamp("2025-06-04")
dia = ancha[ancha["fecha"] == hoy]
print({f"{(v - hoy).days} days": len(g)
       for v, g in dia.groupby("vencimiento")})

{'16 days': 48, '44 days': 45, '72 days': 40}

pivot_table leaves one row per strike with the bid and the ask of the call and of the put. That day there are three expiries in the daily chain, at 16, 44 and 72 days, with 48, 45 and 40 strikes that have both legs.

def hasta_dos_ceros(bid):
    """From K0 outward: drop zero bids; stop on
    finding two in a row."""
    nulo = bid <= 0
    pares = np.flatnonzero(nulo[:-1] & nulo[1:])
    fin = pares[0] if len(pares) else len(bid)
    usar = ~nulo
    usar[fin:] = False
    return usar


def var_cboe(g, T, R, detalle=False):
    K = g["strike"].to_numpy()
    c = ((g["bid_C"] + g["ask_C"]) / 2).to_numpy()
    p = ((g["bid_P"] + g["ask_P"]) / 2).to_numpy()
    j = np.argmin(np.abs(c - p))
    F = K[j] + math.exp(R * T) * (c[j] - p[j])
    i0 = np.flatnonzero(K <= F)[-1]
    ip = np.arange(i0 - 1, -1, -1)       # puts, downward
    ic = np.arange(i0 + 1, len(K))       # calls, upward
    ip = ip[hasta_dos_ceros(g["bid_P"].to_numpy()[ip])][::-1]
    ic = ic[hasta_dos_ceros(g["bid_C"].to_numpy()[ic])]
    idx = np.r_[ip, i0, ic]
    Q = np.r_[p[ip], (c[i0] + p[i0]) / 2, c[ic]]
    Ks = K[idx]
    dK = np.gradient(Ks)   # centered; at the edges, the neighbor
    aporte = 2 / T * dK / Ks ** 2 * Q * math.exp(R * T)
    s2 = aporte.sum() - (F / K[i0] - 1) ** 2 / T
    if detalle:
        return s2, F, K[i0], pd.DataFrame(
            {"strike": Ks, "tipo": ["P"] * len(ip) + ["PC"]
             + ["C"] * len(ic), "Q": Q, "aporte": aporte})
    return s2


for v, g in dia.groupby("vencimiento"):
    n = (v - hoy).days
    T = n / 365
    s2, F, K0, tabla = var_cboe(g, T, tipo(hoy, T), detalle=True)
    puts = tabla.loc[tabla["tipo"] == "P", "aporte"].sum()
    print(f"{n} days: F = {F:.2f}, K0 = {K0:.0f}, "
          f"{len(tabla)} strikes, σ² = {s2:.6f}, "
          f"vol = {100 * math.sqrt(s2):.2f}, "
          f"puts {puts / tabla['aporte'].sum():.0%}")

16 days: F = 5143.89, K0 = 5100, 32 strikes, σ² = 0.065072, vol = 25.51, puts 56%
44 days: F = 5155.03, K0 = 5150, 40 strikes, σ² = 0.067723, vol = 26.02, puts 69%
72 days: F = 5166.23, K0 = 5150, 39 strikes, σ² = 0.069099, vol = 26.29, puts 69%

For each expiry, the forward comes from parity at the strike with the smallest $|C - P|$, $K_0$ is the strike just below and the selection leaves 32, 40 and 39 strikes: in the 16-day one the wings with zero bids are left out. Puts contribute 56%, 69% and 69% of the variance; with skew, the low part of the smile weighs more (and the $1/K^2$ weights favor it even more). The volatility of each expiry, 25.5%, 26.0% and 26.3%, is that of a variance swap at that tenor, not the ATM one: the 44-day ATM of that day was 22.2% (lesson 3.4.5).

def vorx(fecha, grupos):
    """30-day index with the first two expiries of 8 days
    or more (Cboe white paper, monthly version)."""
    vencs = sorted(v for v in grupos if (v - fecha).days >= 8)[:2]
    n1, n2 = [(v - fecha).days for v in vencs]
    s2 = [var_cboe(grupos[v], n / 365, tipo(fecha, n / 365))
          for v, n in zip(vencs, (n1, n2))]
    w = (n1 / 365 * s2[0] * (n2 - 30) / (n2 - n1)
         + n2 / 365 * s2[1] * (30 - n1) / (n2 - n1))
    return 100 * math.sqrt(max(w * 365 / 30, 0.0))


grupos_hoy = dict(tuple(dia.groupby("vencimiento")))
print(f"Computed VORX: {vorx(hoy, grupos_hoy):.2f}")
print(f"Published VORX: {vx.loc[hoy, 'vorx']:.2f}")

Computed VORX: 25.89
Published VORX: 25.89

With $N_1 = 16$ and $N_2 = 44$, the 30 days fall halfway and the two expiries weigh the same in total variance. The result, 25.89, matches the published one.

The figure breaks down the variance of the 44-day expiry by strike. Puts extend down to 2,400 points, 53% below the index, with small but not negligible contributions; calls end much earlier, because the far ones have a zero bid. The $K_0$ bar is the tallest because the put and the call (their average) add up there with the weight of a strike close to the money.

Twelve years recomputed

grupos = {}
for (f, v), g in ancha.groupby(["fecha", "vencimiento"]):
    grupos.setdefault(f, {})[v] = g
calc = pd.Series({f: vorx(f, gv) for f, gv in grupos.items()})
dif = (calc.round(2) - vx["vorx"].reindex(calc.index)).abs()
print(f"{len(calc)} days; maximum difference: {dif.max():.3f}")
sin = vx.index.difference(calc.index)
print("Index days without a chain:", [d.date() for d in sin])

3017 days; maximum difference: 0.000
Index days without a chain: [datetime.date(2019, 7, 23)]

The 3,017 days with a chain give the published VORX to the cent: it is the known answer of the dataset, and the proof that the code applies each rule correctly. The only index day missing, 2019-07-23, is the dia_sin_cadena incident of lesson 3.6.4: that day the vendor did not deliver the chain and the index cannot be recomputed.

ver = pd.read_parquet("opciones3_verdad.parquet")


def atm30(g):
    """True 30-day ATM vol: θ interpolated, θ(0) = 0."""
    th = np.interp(30 / 365, np.r_[0, g["dias"] / 365],
                   np.r_[0, g["theta"]])
    return 100 * math.sqrt(th * 365 / 30)


vx["atm30"] = ver.groupby("fecha").apply(atm30)
d_mod = vx["vorx"] - vx["vorx_modelo"]
d_atm = vx["vorx"] - vx["atm30"]
print(f"VORX − model swap: mean {d_mod.mean():.2f}, "
      f"min {d_mod.min():.2f}, max {d_mod.max():.2f}")
print(f"VORX − 30-day ATM: mean {d_atm.mean():.2f}, "
      f"median {d_atm.median():.2f}")
vx[["vorx", "vorx_modelo", "atm30"]].describe().round(2) \
    .loc[["mean", "50%", "min", "max"]]

VORX − model swap: mean 0.17, min -2.36, max 1.39
VORX − 30-day ATM: mean 2.86, median 2.58
       vorx  vorx_modelo  atm30
mean  19.80        19.63  16.94
50%   17.66        17.37  15.13
min    9.55         9.72   8.50
max   73.98        74.16  63.02

The VORX is on average 0.17 points above the variance swap of the true surface (vorx_modelo), with deviations from −2.36 to 1.39: it is the error of discretizing the integral over strikes, truncating the wings and interpolating between two expiries. Against the 30-day ATM volatility, instead, it is 2.86 points above on average (median 2.58): the difference is not an error, it is the skew. A VIX-like index does not measure ATM volatility, but the square root of an average of the variances of the whole smile, and with skew that average exceeds the at-the-money one.

Realized volatility and risk premium

The fair comparison is today's implied against the realized afterward: the VORX of day $t$ looks at the following 30 calendar days, about 21 sessions. Using the realized of the previous 21 sessions would be comparing today's price with a figure that was already known.

r = np.log(ind["cierre"]).diff()
real_post = 100 * np.sqrt(252 * (r ** 2).rolling(21).mean()) \
    .shift(-21)                       # sessions t+1 … t+21
real_prev = 100 * np.sqrt(252 * (r ** 2).rolling(21).mean())
prima = (vx["vorx"] - real_post).dropna()
print(f"VORX − realized premium: mean {prima.mean():.2f}, "
      f"median {prima.median():.2f}")
print(f"Days with realized above: {(prima < 0).mean():.1%}")
print(f"Worst day: {prima.idxmin().date()}, {prima.min():.1f} points")
p_atm = (vx["atm30"] - real_post).dropna()
print(f"With the 30-day ATM: mean {p_atm.mean():.2f}, "
      f"median {p_atm.median():.2f}")
var_prima = (vx["vorx"] ** 2 - real_post ** 2).dropna() / 100
print(f"In variance (points² / 100): mean {var_prima.mean():.2f}")

VORX − realized premium: mean 5.30, median 6.38
Days with realized above: 11.1%
Worst day: 2020-03-06, -76.5 points
With the 30-day ATM: mean 2.43, median 3.88
In variance (points² / 100): mean 1.53

The mean premium is 5.30 volatility points and the median, 6.38. Realized exceeds the VORX on 11.1% of the days, and the worst case is March 6, 2020: the VORX stood at 23.2 and the realized of the following 21 sessions, which include the double-digit falls of March 12 and 16, was 99.7%, a difference of −76.5 points. With the 30-day ATM volatility, the premium is smaller (2.43 on average) because part of what the VORX charges is the skew. In variance, which is what is traded in a swap, the mean difference is 153 squared volatility points (0.0153 of annual variance).

In the data the premium is planted: the implied variance is between 1.0 and 1.3 times the expected one (more in calm than in stress), according to the data README. The academic reference on the variance premium in real markets is Carr and Wu (2009), in the resources; no real figure is given here.

anual = prima.groupby(prima.index.year).agg(
    media="mean", mediana="median", minimo="min",
    negativa=lambda z: (z < 0).mean())
anual.round(2).T

fecha     2014   2015   2016  2017   2018  2019   2020  2021  2022  2023  2024   2025
media     3.91   4.38   4.19  4.26   6.04  6.90   5.08  7.18  5.67  6.35  5.41   4.18
mediana   4.23   5.68   5.50  5.02   8.13  7.07  11.92  8.62  6.37  7.15  5.63   6.33
minimo   -4.19 -15.67 -16.19 -8.09 -16.95 -1.40 -76.46 -7.03 -9.44 -1.15 -2.88 -24.94
negativa  0.11   0.22   0.15  0.17   0.16  0.02   0.10  0.09  0.08  0.04  0.06   0.14

The median is positive every year, from 4.2 in 2014 to 11.9 in 2020; the fraction of negative days ranges from 2% (2019) to 22% (2015). The minimums tell another story: −76.5 in 2020, −24.9 in 2025, −17.0 in 2018 and −16.2 in 2016, the years with stress periods in the data script. The premium is stable in the center and very unstable in the tail.

The figure shows the whole distribution, with the number of days on a log scale so that the tail is visible: the mass between 0 and 15 points and, to the left of zero (highlighted), a long tail with an isolated group between −76 and −62 that is the 16 sessions from February 13 to March 6, 2020.

import statsmodels.api as sm

d = pd.concat({"real_post": real_post, "vorx": vx["vorx"],
               "real_prev": real_prev}, axis=1).dropna()
d = d.iloc[::21]          # no overlap: one every 21 sessions
for x_cols in (["vorx"], ["real_prev"], ["vorx", "real_prev"]):
    m = sm.OLS(d["real_post"], sm.add_constant(d[x_cols])).fit()
    coef = ", ".join(f"{c} {m.params[c]:.2f} (t {m.tvalues[c]:.1f})"
                     for c in x_cols)
    print(f"{coef} · R² {m.rsquared:.2f} · n {int(m.nobs)}")

vorx 0.70 (t 7.9) · R² 0.31 · n 142
real_prev 0.41 (t 5.3) · R² 0.17 · n 142
vorx 0.91 (t 5.5), real_prev -0.19 (t -1.5) · R² 0.32 · n 142

With one observation every 21 sessions (142, with no overlap), the VORX explains 31% of the variance of the subsequent realized, against 17% for the past realized; together, 32%, and the past one loses its significance (t = −1.5). The VORX contains the information of past volatility and something more. Its slope of 0.70 (alone) says that it is an upward-biased predictor: it carries the premium inside.

The VSTOXX and the European index of the data

The VSTOXX is the volatility index of the EURO STOXX 50, developed by Goldman Sachs and Deutsche Börse and calculated by STOXX with Eurex OESX options. Its variance per expiry uses the same formula as the VIX, but the STOXX guide22Note 2 · STOXX, Strategy Index Guide, September 2024 edition, section 9 (EURO STOXX 50 Volatility). Link in the resources. changes quite a few details:

The practical consequence is that the two indices are not interchangeable even though they move together: with wide spreads (stress), the VSTOXX can use settlement or last-trade prices where the VIX uses mids, and the 0.5-point filter trims its wings earlier. In the course data, the fictitious European index VGX has its volatility index VVGX, computed as the square root of the 30-day variance swap of its own model (there is no published VGX chain, so it cannot be recomputed):

rg = np.log(vx["vgx"]).diff()
real_post_g = 100 * np.sqrt(252 * (rg ** 2).rolling(21).mean()) \
    .shift(-21)
prima_g = (vx["vvgx"] - real_post_g).dropna()
print(f"VVGX: mean {vx['vvgx'].mean():.1f}, "
      f"VORX: mean {vx['vorx'].mean():.1f}")
print(f"VVGX premium: mean {prima_g.mean():.2f}, median "
      f"{prima_g.median():.2f}, negative {(prima_g < 0).mean():.1%}")
print(f"ORX-VGX return correlation: {r.corr(rg):.2f}")
print("VORX-VVGX daily change correlation: "
      f"{vx['vorx'].diff().corr(vx['vvgx'].diff()):.2f}")
print("Correlation of ORX with VORX changes: "
      f"{r.corr(vx['vorx'].diff()):.2f}")

VVGX: mean 21.3, VORX: mean 19.8
VVGX premium: mean 5.41, median 6.16, negative 11.0%
ORX-VGX return correlation: 0.82
VORX-VVGX daily change correlation: 0.96
Correlation of ORX with VORX changes: -0.61

The VVGX has a mean of 21.3 against 19.8 for the VORX, a similar premium (5.41 mean, 6.16 median, negative on 11.0% of the days) and daily changes almost identical to those of the VORX (correlation 0.96), because in the data model its variance is 1.2 times that of ORX with a small noise of its own. In real markets the relationship between the VIX and the VSTOXX is close, but not that close: each has its own options market, its trading hours and its rules. The negative correlation of −0.61 between the ORX return and the VORX change is the leverage effect: implied volatility rises when the index falls.

Exercises

1. Compute the VORX for Wednesday, June 18, 2025, when the first expiry is two days away. Which expiries enter and why? Check the result against the published one.

f = pd.Timestamp("2025-06-18")
print({(v - f).days: len(g) for v, g in grupos[f].items()})
print(f"Computed {vorx(f, grupos[f]):.2f}, "
      f"published {vx.loc[f, 'vorx']:.2f}")

{2: 48, 30: 45, 58: 40}
Computed 25.11, published 25.11

2. Repeat the premium with the realized of the previous 21 sessions instead of the following ones. How much do the mean and the fraction of negative days change? Why is that version useless for evaluating a strategy?

prima_prev = (vx["vorx"] - real_prev).dropna()
print(f"Mean {prima_prev.mean():.2f} (subsequent {prima.mean():.2f})")
print(f"Negative {(prima_prev < 0).mean():.1%} "
      f"(subsequent {(prima < 0).mean():.1%})")
print(f"Minimum {prima_prev.min():.1f} on {prima_prev.idxmin().date()}")

Mean 5.33 (subsequent 5.30)
Negative 12.3% (subsequent 11.1%)
Minimum -36.6 on 2020-04-06

3. In the chain of June 4, 2025, the VORX uses 32 of the 48 strikes of the 16-day expiry. What would happen to its variance if all of them were used, with the mid of the zero-bid options? Reason it out with the formula and check it.

g = grupos_hoy[pd.Timestamp("2025-06-20")]
T = 16 / 365
R = tipo(hoy, T)
K = g["strike"].to_numpy()
c = ((g["bid_C"] + g["ask_C"]) / 2).to_numpy()
p = ((g["bid_P"] + g["ask_P"]) / 2).to_numpy()
j = np.argmin(np.abs(c - p))
F = K[j] + math.exp(R * T) * (c[j] - p[j])
i0 = np.flatnonzero(K <= F)[-1]
Q = np.r_[p[:i0], (c[i0] + p[i0]) / 2, c[i0 + 1:]]
todos = (2 / T * np.gradient(K) / K ** 2 * Q
         * math.exp(R * T)).sum() - (F / K[i0] - 1) ** 2 / T
regla = var_cboe(g, T, R)
print(f"Cboe rule: σ² = {regla:.6f} "
      f"(vol {100 * math.sqrt(regla):.2f})")
print(f"All strikes: σ² = {todos:.6f} "
      f"(vol {100 * math.sqrt(todos):.2f})")

Cboe rule: σ² = 0.065072 (vol 25.51)
All strikes: σ² = 0.065269 (vol 25.55)