Measuring the Variance Risk Premium

Published on 09/2026

Is This Option Expensive?

If you're looking at one option contract, the obvious question is whether it's priced richly or cheaply. The naive way to answer that is to compare its implied volatility (the volatility the market is charging for) to how volatile you think the stock will actually be, and call the difference your profit margin.

The problem is, that answer is wrong. Option prices sit above the volatility that actually turns up almost all the time - selling an option is basically selling insurance, and insurance is normally sold above cost. So finding a gap between implied and realised volatility isn't a discovery, it's the default state. The real question is whether today's gap is bigger than this particular stock's usual gap, and that's what this tool works out.

Now, decomposing the problem:

Getting the data
Measuring realised volatility
Building the historical premium
Putting an honest error bar on it
Forecasting volatility
Cleaning up today's price

Removing the earnings bump
Accounting for the strike

Scoring the trade
Testing it
Limitations

Getting the Data

I used the Interactive Brokers API (through the ib_async library) for the data. The first idea most people have is to pull the price history of the option contract itself, but that doesn't work: a fixed strike slides around as the share price moves. Nvidia's $160 put was far out-of-the-money in February (market price ~$175-$195 throughout the month) and much closer to the money at the end of March (market price ~$163-$173 in the final days), so averaging that history describes nothing in particular.

Instead, we pull two histories for the stock: five years of daily open, high, low and close prices, and five years of the stock's typical at-the-money implied volatility, which IBKR provides as a single series. That's only two requests, which matters because IBKR only allows about 60 historical requests every 10 minutes. The specific strike we care about then gets treated as an adjustment on top, later on.

python

def daily_bars(self, symbol, duration="5 Y"):
    bars = self.ib.reqHistoricalData(
        self._qualified_stock(symbol), endDateTime="", durationStr=duration,
        barSizeSetting="1 day", whatToShow="TRADES", useRTH=True, formatDate=1)
    df = pd.DataFrame([b.__dict__ for b in bars])
    df["date"] = pd.to_datetime(df["date"])
    return df.set_index("date")[["open", "high", "low", "close", "volume"]]

def iv_history(self, symbol, duration="5 Y"):
    bars = self.ib.reqHistoricalData(
        self._qualified_stock(symbol), endDateTime="", durationStr=duration,
        barSizeSetting="1 day", whatToShow="OPTION_IMPLIED_VOLATILITY",
        useRTH=True, formatDate=1)
    df = pd.DataFrame([b.__dict__ for b in bars])
    df["date"] = pd.to_datetime(df["date"])
    return df.set_index("date")["close"].rename("iv")

python

def daily_bars(self, symbol, duration="5 Y"):
    bars = self.ib.reqHistoricalData(
        self._qualified_stock(symbol), endDateTime="", durationStr=duration,
        barSizeSetting="1 day", whatToShow="TRADES", useRTH=True, formatDate=1)
    df = pd.DataFrame([b.__dict__ for b in bars])
    df["date"] = pd.to_datetime(df["date"])
    return df.set_index("date")[["open", "high", "low", "close", "volume"]]

def iv_history(self, symbol, duration="5 Y"):
    bars = self.ib.reqHistoricalData(
        self._qualified_stock(symbol), endDateTime="", durationStr=duration,
        barSizeSetting="1 day", whatToShow="OPTION_IMPLIED_VOLATILITY",
        useRTH=True, formatDate=1)
    df = pd.DataFrame([b.__dict__ for b in bars])
    df["date"] = pd.to_datetime(df["date"])
    return df.set_index("date")["close"].rename("iv")

Measuring Realised Volatility

The obvious way to measure how volatile a stock has been is to take the closing prices and see how much they move from day to day. However, that throws away most of the information - we also have the open, high and low for each day. It also can't tell apart a stock that gaps overnight (closes at 100 and opens at 94) from one that drifts down slowly, and single stocks gap on news all the time.

The Yang-Zhang estimator fixes this by splitting each day into three pieces: the overnight gap, the move from open to close, and a measure built from the high and low (the Rogers-Satchell term). For day $t$ , with open $O_t$ , high $H_t$ , low $L_t$ and close $C_t$ , those are:

o_t = \ln\frac{O_t}{C_{t-1}}, \qquad c_t = \ln\frac{C_t}{O_t}, \qquad rs_t = \ln\frac{H_t}{C_t}\ln\frac{H_t}{O_t} + \ln\frac{L_t}{C_t}\ln\frac{L_t}{O_t}

Over a window of $n$ days, we take the sample variance of the overnight and open-to-close moves, and the average of the Rogers-Satchell term, then weight them together:

\sigma^2_{YZ} = \sigma^2_{o} + k\,\sigma^2_{c} + (1-k)\,\sigma^2_{RS}, \qquad k = \frac{0.34}{1.34 + \frac{n+1}{n-1}}

The coefficient $k$ is chosen to make the estimate as stable as possible, and the result is annualised by multiplying by the number of trading days in a year and taking the square root. This is roughly 5-14x more efficient than the closing-price method on the same data, so we can use a shorter window and still get a reliable number.

python

def yang_zhang(df, window):
    o, h, l, c = df["open"], df["high"], df["low"], df["close"]

    log_oc_prev = np.log(o / c.shift(1))      # overnight gap
    log_co = np.log(c / o)                    # intraday move
    rs = np.log(h / c) * np.log(h / o) + np.log(l / c) * np.log(l / o)

    var_overnight = log_oc_prev.rolling(window).var(ddof=1)
    var_open_close = log_co.rolling(window).var(ddof=1)
    var_rs = rs.rolling(window).mean()

    k = 0.34 / (1.34 + (window + 1) / (window - 1))
    var_yz = var_overnight + k * var_open_close + (1 - k) * var_rs
    return np.sqrt(var_yz.clip(lower=0) * TRADING_DAYS)

python

def yang_zhang(df, window):
    o, h, l, c = df["open"], df["high"], df["low"], df["close"]

    log_oc_prev = np.log(o / c.shift(1))      # overnight gap
    log_co = np.log(c / o)                    # intraday move
    rs = np.log(h / c) * np.log(h / o) + np.log(l / c) * np.log(l / o)

    var_overnight = log_oc_prev.rolling(window).var(ddof=1)
    var_open_close = log_co.rolling(window).var(ddof=1)
    var_rs = rs.rolling(window).mean()

    k = 0.34 / (1.34 + (window + 1) / (window - 1))
    var_yz = var_overnight + k * var_open_close + (1 - k) * var_rs
    return np.sqrt(var_yz.clip(lower=0) * TRADING_DAYS)

This estimator gets used twice - once to build the history of what volatility actually turned out to be, and again as the input to the forecasting model - so a mistake here would affect every number the tool prints.

Building the Historical Premium

For every past day, we take what implied volatility was on that day and pair it with how volatile the stock actually turned out to be over the following month (21 trading days). The difference is that day's insurance premium - the variance risk premium. Writing $h = 21$ for the horizon:

VRP_t = IV_t - RV_{t \to t+h}

where $RV_{t \to t+h}$ is the Yang-Zhang volatility over days $t+1$ to $t+h$ .

The forward direction is the most important detail in the whole project. If we compared implied volatility to the previous month instead, we'd be measuring something completely different - that calm periods follow calm periods and chaotic ones follow chaotic ones - and it would look convincing while being meaningless. Shifting the trailing estimate back by the horizon puts the future window on each row:

python

def forward_realised_vol(df, horizon, estimator=yang_zhang):
    trailing = estimator(df, horizon)
    # trailing.loc[t] covers (t-horizon, t]; shifting back by `horizon` puts the
    # window (t, t+horizon] on the row indexed t.
    return trailing.shift(-horizon)

def build_vrp_series(iv_history, ohlc, horizon=21):
    rv_fwd = forward_realised_vol(ohlc, horizon)
    out = pd.DataFrame({"iv": iv_history, "rv_forward": rv_fwd}).dropna()
    out["vrp"] = out["iv"] - out["rv_forward"]
    return out

python

def forward_realised_vol(df, horizon, estimator=yang_zhang):
    trailing = estimator(df, horizon)
    # trailing.loc[t] covers (t-horizon, t]; shifting back by `horizon` puts the
    # window (t, t+horizon] on the row indexed t.
    return trailing.shift(-horizon)

def build_vrp_series(iv_history, ohlc, horizon=21):
    rv_fwd = forward_realised_vol(ohlc, horizon)
    out = pd.DataFrame({"iv": iv_history, "rv_forward": rv_fwd}).dropna()
    out["vrp"] = out["iv"] - out["rv_forward"]
    return out

Putting an Honest Error Bar on It

Averaging those daily premiums over roughly three years (756 trading days) gives the stock's normal premium. That sounds like a lot of evidence, but it isn't: Monday's reading looks at the next 21 days, and Tuesday's reading looks at 20 of the same days plus one new one. The two are almost the same observation, so there are far fewer independent readings than there are rows.

If we calculated the standard error the normal way, it would come out roughly 3-5 times too small, and a "normal premium" that's really just noise would look statistically solid. The Newey-West correction fixes this by looking at how strongly each reading correlates with its neighbours (1 day apart, 2 days apart, all the way to 21) and adding those overlaps back into the uncertainty, with closer neighbours weighted more heavily.

Written out, with $n$ premium readings $x_1, \dots, x_n$ , their mean $\bar{x}$ and deviations $e_t = x_t - \bar{x}$ , the covariance at each lag $L$ and its weight are:

\hat{\gamma}_L = \frac{1}{n} \sum_{t=L+1}^{n} e_t\, e_{t-L}, \qquad w_L = 1 - \frac{L}{h+1}

and the standard error of the mean becomes:

\mathrm{SE}_{NW}(\bar{x}) = \sqrt{\frac{1}{n} \left( \hat{\gamma}_0 + 2 \sum_{L=1}^{h} w_L\, \hat{\gamma}_L \right)}

The ordinary standard error is just the $\hat{\gamma}_0$ term on its own - the sum is the correction for the overlapping windows.

python

def newey_west_se(x, lags):
    n = len(x)
    e = x - x.mean()
    s = float(e @ e) / n
    for L in range(1, min(lags, n - 1) + 1):
        w = 1.0 - L / (lags + 1.0)
        s += 2.0 * w * float(e[L:] @ e[:-L]) / n
    return float(np.sqrt(max(s, 0.0) / n))

def vrp_baseline(vrp_df, lookback_days=756, horizon=21):
    s = vrp_df["vrp"]
    if lookback_days:
        s = s.iloc[-lookback_days:]
    v = s.dropna().values
    return {
        "n_obs": len(v),
        "mean": float(np.mean(v)),
        "median": float(np.median(v)),
        "std": float(np.std(v, ddof=1)),
        "se_newey_west": newey_west_se(v, lags=horizon),
        "pct_positive": float(np.mean(v > 0)),
        "p10": float(np.percentile(v, 10)),
        "p90": float(np.percentile(v, 90)),
    }

python

def newey_west_se(x, lags):
    n = len(x)
    e = x - x.mean()
    s = float(e @ e) / n
    for L in range(1, min(lags, n - 1) + 1):
        w = 1.0 - L / (lags + 1.0)
        s += 2.0 * w * float(e[L:] @ e[:-L]) / n
    return float(np.sqrt(max(s, 0.0) / n))

def vrp_baseline(vrp_df, lookback_days=756, horizon=21):
    s = vrp_df["vrp"]
    if lookback_days:
        s = s.iloc[-lookback_days:]
    v = s.dropna().values
    return {
        "n_obs": len(v),
        "mean": float(np.mean(v)),
        "median": float(np.median(v)),
        "std": float(np.std(v, ddof=1)),
        "se_newey_west": newey_west_se(v, lags=horizon),
        "pct_positive": float(np.mean(v > 0)),
        "p10": float(np.percentile(v, 10)),
        "p90": float(np.percentile(v, 90)),
    }

Forecasting Volatility

Next we need a forecast of how volatile the stock will actually be over the option's life. For this I used HAR-RV (the Heterogeneous Autoregressive model of Realised Volatility, from Corsi in 2009). The idea is that different traders work on different timescales - day traders, people rebalancing weekly, and institutions thinking in months - so the model feeds in volatility over three horizons at once: the last day, the last 5 days and the last 22 days.

It's just a straight-line regression, fitted in logarithms because volatility is heavily skewed (it can spike to four times normal but can never go below zero). Despite being this simple, it's genuinely hard to beat with more complicated models. I also calculate an exponentially-weighted average of the volatility, but only as a sanity check - if the two disagree wildly, something is probably wrong with the data.

With $RV_t$ the trailing Yang-Zhang volatility on day $t$ , the weekly and monthly inputs are just its averages over the last 5 and 22 days:

RV^{(w)}_t = \frac{1}{5} \sum_{i=0}^{4} RV_{t-i}, \qquad RV^{(m)}_t = \frac{1}{22} \sum_{i=0}^{21} RV_{t-i}

and the regression predicts the (log) volatility over the next month from all three:

\ln RV_{t \to t+h} = \beta_0 + \beta_d \ln RV_t + \beta_w \ln RV^{(w)}_t + \beta_m \ln RV^{(m)}_t + \varepsilon_t

The $\beta$ s are fitted by least squares, and the forecast is today's fitted value taken back out of logs:

\hat{\sigma} = \exp\left( \hat{\beta}_0 + \hat{\beta}_d \ln RV_t + \hat{\beta}_w \ln RV^{(w)}_t + \hat{\beta}_m \ln RV^{(m)}_t \right)

python

class HARRV:
    def __init__(self, horizon=21):
        self.horizon = horizon
        self.coef_ = None
        self.resid_std_ = None

    @staticmethod
    def _features(rv_daily):
        return pd.DataFrame({
            "rv_d": np.log(rv_daily),
            "rv_w": np.log(rv_daily.rolling(5).mean()),
            "rv_m": np.log(rv_daily.rolling(22).mean()),
        })

    def fit(self, rv_daily, rv_forward):
        X = self._features(rv_daily)
        y = np.log(rv_forward)
        data = pd.concat([X, y.rename("y")], axis=1).replace(
            [np.inf, -np.inf], np.nan).dropna()
        if len(data) < 60:
            raise ValueError(f"only {len(data)} usable rows; need >= 60")

        A = np.column_stack([np.ones(len(data)), data[["rv_d", "rv_w", "rv_m"]].values])
        self.coef_, *_ = np.linalg.lstsq(A, data["y"].values, rcond=None)
        self.resid_std_ = float(np.std(data["y"].values - A @ self.coef_, ddof=4))
        return self

    def predict(self, rv_daily):
        if self.coef_ is None:
            raise RuntimeError("call fit() first")
        x = self._features(rv_daily).iloc[-1]
        if x.isna().any():
            raise ValueError("insufficient history for features")
        return float(np.exp(self.coef_ @ np.r_[1.0, x.values]))

python

class HARRV:
    def __init__(self, horizon=21):
        self.horizon = horizon
        self.coef_ = None
        self.resid_std_ = None

    @staticmethod
    def _features(rv_daily):
        return pd.DataFrame({
            "rv_d": np.log(rv_daily),
            "rv_w": np.log(rv_daily.rolling(5).mean()),
            "rv_m": np.log(rv_daily.rolling(22).mean()),
        })

    def fit(self, rv_daily, rv_forward):
        X = self._features(rv_daily)
        y = np.log(rv_forward)
        data = pd.concat([X, y.rename("y")], axis=1).replace(
            [np.inf, -np.inf], np.nan).dropna()
        if len(data) < 60:
            raise ValueError(f"only {len(data)} usable rows; need >= 60")

        A = np.column_stack([np.ones(len(data)), data[["rv_d", "rv_w", "rv_m"]].values])
        self.coef_, *_ = np.linalg.lstsq(A, data["y"].values, rcond=None)
        self.resid_std_ = float(np.std(data["y"].values - A @ self.coef_, ddof=4))
        return self

    def predict(self, rv_daily):
        if self.coef_ is None:
            raise RuntimeError("call fit() first")
        x = self._features(rv_daily).iloc[-1]
        if x.isna().any():
            raise ValueError("insufficient history for features")
        return float(np.exp(self.coef_ @ np.r_[1.0, x.values]))

Cleaning Up Today's Price

Before comparing today's implied volatility against the historical premium, it needs to be the same kind of number. There are two adjustments for this.

Removing the Earnings Bump

An option that expires shortly after an earnings date carries the price of that one event, so it would look hugely overpriced when really it's just sitting before an earnings date.

To fix this, we quote a second, longer-dated option too. Total variance over an option's life (with $T$ in years) is the normal day-to-day variance plus the one-off earnings jump $J$ , and dividing through by $T$ gives the implied volatility we actually see:

\sigma_{total}^2\, T = \sigma_{base}^2\, T + J^2 \quad \Longrightarrow \quad \sigma_{total}^2 = \sigma_{base}^2 + \frac{J^2}{T}

The jump gets diluted more in the longer option (it's divided by a bigger $T$ ), so comparing the two expiries lets us work out how big a move the market expects on earnings and subtract it. Writing that equation for the front ( $f$ ) and back ( $b$ ) expiries and subtracting one from the other cancels $\sigma_{base}$ out:

\sigma_f^2 - \sigma_b^2 = J^2 \left( \frac{1}{T_f} - \frac{1}{T_b} \right) \quad \Longrightarrow \quad J = \sqrt{\frac{\sigma_f^2 - \sigma_b^2}{1/T_f - 1/T_b}}

Then putting $J$ back into the front-month equation strips the earnings out of today's price:

\sigma_{base} = \sqrt{\sigma_f^2 - \frac{J^2}{T_f}}

If the longer-dated option has the higher implied volatility, the top of that fraction for $J$ is negative - there's no earnings hump to remove, so the function refuses rather than inventing one.

python

def implied_jump_from_term_structure(iv_front, T_front, iv_back, T_back):
    num = iv_front**2 - iv_back**2
    den = (1.0 / T_front) - (1.0 / T_back)
    if den <= 0 or num <= 0:
        raise ValueError("term structure not in event-hump configuration")
    return float(np.sqrt(num / den))

def strip_event_variance(iv_total, T_years, implied_jump):
    base_var = iv_total**2 - (implied_jump**2) / T_years
    if base_var <= 0:
        raise ValueError("implied jump exceeds total variance; check inputs")
    return float(np.sqrt(base_var))

python

def implied_jump_from_term_structure(iv_front, T_front, iv_back, T_back):
    num = iv_front**2 - iv_back**2
    den = (1.0 / T_front) - (1.0 / T_back)
    if den <= 0 or num <= 0:
        raise ValueError("term structure not in event-hump configuration")
    return float(np.sqrt(num / den))

def strip_event_variance(iv_total, T_years, implied_jump):
    base_var = iv_total**2 - (implied_jump**2) / T_years
    if base_var <= 0:
        raise ValueError("implied jump exceeds total variance; check inputs")
    return float(np.sqrt(base_var))

Accounting for the Strike

The historical premium was built at the money, but an out-of-the-money put carries an extra premium just for being out of the money (the volatility skew). If we ignored this, every out-of-the-money put would look mispriced. So we take the difference between our option's implied volatility and the at-the-money one, and add it to the expected side of the comparison rather than counting it as a finding. For an option with strike $K$ , that spread is:

s = IV_K - IV_{ATM}

Scoring the Trade

Now we can put it all together. The tool prints three numbers:

Gap - today's implied volatility minus the forecast. This is what a basic screener would show you.
Typical gap - what that gap normally is for this stock, plus the strike adjustment.
Edge - the gap minus the typical gap. This is the only number that actually carries information.

In symbols, with $\hat{\sigma}$ the HAR-RV forecast, $\overline{VRP}$ the stock's average historical premium and $s$ the strike spread:

\text{gap} = IV_K - \hat{\sigma}, \qquad \text{typical} = \overline{VRP} + s, \qquad \text{edge} = \text{gap} - \text{typical}

If the edge is smaller than twice the Newey-West error bar, the verdict is "no edge", since we're inside the noise of the baseline itself. The z-score printed next to it measures the edge against how much the premium normally swings around (its standard deviation), rather than against the error bar:

|\text{edge}| < 2\, \mathrm{SE}_{NW} \;\Longrightarrow\; \text{no edge}, \qquad z = \frac{\text{edge}}{\operatorname{sd}(VRP)}

python

def score_trade(market_iv, rv_forecast, baseline, strike_skew_spread=0.0):
    gap = market_iv - rv_forecast
    typical = baseline["mean"] + strike_skew_spread
    edge = gap - typical
    z = edge / baseline["std"] if baseline["std"] > 0 else float("nan")

    if abs(edge) < 2 * baseline["se_newey_west"]:
        verdict = "NO EDGE — inside the noise band of the baseline itself"
    elif edge > 0:
        verdict = "Overpaid to sell vol relative to the historical premium"
    else:
        verdict = "Underpaid to sell / relatively cheap to buy"

    return {
        "market_iv": market_iv,
        "rv_forecast": rv_forecast,
        "gap": gap,
        "typical_gap": typical,
        "edge": edge,
        "z_score": z,
        "baseline_se": baseline["se_newey_west"],
        "verdict": verdict,
    }

python

def score_trade(market_iv, rv_forecast, baseline, strike_skew_spread=0.0):
    gap = market_iv - rv_forecast
    typical = baseline["mean"] + strike_skew_spread
    edge = gap - typical
    z = edge / baseline["std"] if baseline["std"] > 0 else float("nan")

    if abs(edge) < 2 * baseline["se_newey_west"]:
        verdict = "NO EDGE — inside the noise band of the baseline itself"
    elif edge > 0:
        verdict = "Overpaid to sell vol relative to the historical premium"
    else:
        verdict = "Underpaid to sell / relatively cheap to buy"

    return {
        "market_iv": market_iv,
        "rv_forecast": rv_forecast,
        "gap": gap,
        "typical_gap": typical,
        "edge": edge,
        "z_score": z,
        "baseline_se": baseline["se_newey_west"],
        "verdict": verdict,
    }

To see why this matters, imagine two candidates (these numbers are made up just to show the idea). An Nvidia put has a gap of 8 points, but Nvidia's typical gap is 5.5, so its edge is +2.5. A Coca-Cola put has a smaller gap of 6.2, but its typical gap is only 3.1, so its edge is +3.1.

Sorted by gap, Nvidia wins comfortably. Sorted by edge, Coca-Cola is the better trade - Nvidia only looks dramatic because it's a volatile stock, where an 8-point gap is close to ordinary. Ranking by the raw gap steers you towards volatile stocks that are priced exactly as they should be.

The last function wires everything together. I kept all the statistics separate from the IBKR code on purpose, so this part can run and be tested without a broker connection:

python

def analyse(iv_history, ohlc, market_iv, atm_iv=None, horizon=HORIZON):
    iv_history, ohlc = iv_history.align(ohlc, join="inner", axis=0)

    # right-hand column: what the premium normally is, at the money
    vrp_df = build_vrp_series(iv_history, ohlc, horizon)
    base = vrp_baseline(vrp_df, lookback_days=756, horizon=horizon)

    # middle column: forecast realised vol over the option's life
    rv_daily = yang_zhang(ohlc, window=horizon)
    rv_fwd = forward_realised_vol(ohlc, horizon)
    model = HARRV(horizon).fit(rv_daily, rv_fwd)
    forecast = model.predict(rv_daily)

    # carry the skew across, since the baseline is an ATM object
    skew_spread = (market_iv - atm_iv) if atm_iv is not None else 0.0

    row = score_trade(market_iv, forecast, base, skew_spread)
    row["_baseline"] = base
    row["_ewma_check"] = ewma_vol(ohlc.tail(120))
    row["_skew_spread"] = skew_spread
    return row

python

def analyse(iv_history, ohlc, market_iv, atm_iv=None, horizon=HORIZON):
    iv_history, ohlc = iv_history.align(ohlc, join="inner", axis=0)

    # right-hand column: what the premium normally is, at the money
    vrp_df = build_vrp_series(iv_history, ohlc, horizon)
    base = vrp_baseline(vrp_df, lookback_days=756, horizon=horizon)

    # middle column: forecast realised vol over the option's life
    rv_daily = yang_zhang(ohlc, window=horizon)
    rv_fwd = forward_realised_vol(ohlc, horizon)
    model = HARRV(horizon).fit(rv_daily, rv_fwd)
    forecast = model.predict(rv_daily)

    # carry the skew across, since the baseline is an ATM object
    skew_spread = (market_iv - atm_iv) if atm_iv is not None else 0.0

    row = score_trade(market_iv, forecast, base, skew_spread)
    row["_baseline"] = base
    row["_ewma_check"] = ewma_vol(ohlc.tail(120))
    row["_skew_spread"] = skew_spread
    return row

Testing It

Since this can run offline, I wrote a self-test that simulates a stock with a known premium of exactly 4 volatility points, runs it through the whole pipeline and checks that it comes back out. It also plants a known earnings jump (6%) and checks the decomposition recovers it exactly.

Simulating a Stock

First we need a fake stock whose true volatility we know. Real volatility wanders around but tends to drift back towards a normal level, so we simulate it as a mean-reverting series around 35%: each day it moves 97% of the way back from wherever it was, plus a small random shock, and it's clipped between 12% and 120% so it never does anything absurd. Using a fixed random seed (7) means the test gives exactly the same numbers every time:

\sigma_t = \operatorname{clip}\big( 0.35 + 0.97\,(\sigma_{t-1} - 0.35) + \eta_t,\; 0.12,\; 1.20 \big), \qquad \eta_t \sim N(0,\, 0.012^2)

python

rng = np.random.default_rng(7)
n = 1500

# mean-reverting true vol around 35%
true_vol = np.zeros(n)
true_vol[0] = 0.35
for i in range(1, n):
    true_vol[i] = np.clip(
        0.35 + 0.97 * (true_vol[i - 1] - 0.35) + rng.normal(0, 0.012),
        0.12, 1.20)

python

rng = np.random.default_rng(7)
n = 1500

# mean-reverting true vol around 35%
true_vol = np.zeros(n)
true_vol[0] = 0.35
for i in range(1, n):
    true_vol[i] = np.clip(
        0.35 + 0.97 * (true_vol[i - 1] - 0.35) + rng.normal(0, 0.012),
        0.12, 1.20)

Then we turn that volatility into daily prices. Yang-Zhang needs all four prices for each day, so we can't just simulate closing prices - we build each day in three pieces, which mirror the three pieces the estimator splits a day into: an overnight gap from yesterday's close to today's open, an intraday move from the open to the close, and a random range that pushes the high above and the low below both of them. Each piece is scaled by that day's true volatility, so a volatile day really does produce bigger moves. We simulate 1,500 business days, which is about 6 years and plenty for a 3-year baseline. Writing $\Delta t = 1/252$ for one trading day, the three random pieces are:

g_t \sim N\big(0,\, (0.4\, \sigma_t \sqrt{\Delta t})^2\big), \qquad u_t \sim N\big(0,\, (0.9\, \sigma_t \sqrt{\Delta t})^2\big), \qquad r_t \sim N\big(0,\, (0.7\, \sigma_t \sqrt{\Delta t})^2\big)

and each day is built from them:

O_t = C_{t-1}\, e^{g_t}, \qquad C_t = O_t\, e^{u_t}, \qquad H_t = \max(O_t, C_t)\, e^{|r_t|}, \qquad L_t = \min(O_t, C_t)\, e^{-|r_t|}

python

dt = 1 / 252
close = np.zeros(n)
close[0] = 100.0
o = np.zeros(n); h = np.zeros(n); lo = np.zeros(n)
for i in range(1, n):
    gap = rng.normal(0, true_vol[i] * np.sqrt(dt) * 0.4)
    o[i] = close[i - 1] * np.exp(gap)
    intraday = rng.normal(0, true_vol[i] * np.sqrt(dt) * 0.9)
    close[i] = o[i] * np.exp(intraday)
    rng_width = abs(rng.normal(0, true_vol[i] * np.sqrt(dt) * 0.7))
    h[i] = max(o[i], close[i]) * np.exp(rng_width)
    lo[i] = min(o[i], close[i]) * np.exp(-rng_width)
o[0], h[0], lo[0] = 100.0, 100.5, 99.5

idx = pd.bdate_range("2020-01-01", periods=n)
ohlc = pd.DataFrame({"open": o, "high": h, "low": lo, "close": close,
                     "volume": 1e6}, index=idx).iloc[1:]

python

dt = 1 / 252
close = np.zeros(n)
close[0] = 100.0
o = np.zeros(n); h = np.zeros(n); lo = np.zeros(n)
for i in range(1, n):
    gap = rng.normal(0, true_vol[i] * np.sqrt(dt) * 0.4)
    o[i] = close[i - 1] * np.exp(gap)
    intraday = rng.normal(0, true_vol[i] * np.sqrt(dt) * 0.9)
    close[i] = o[i] * np.exp(intraday)
    rng_width = abs(rng.normal(0, true_vol[i] * np.sqrt(dt) * 0.7))
    h[i] = max(o[i], close[i]) * np.exp(rng_width)
    lo[i] = min(o[i], close[i]) * np.exp(-rng_width)
o[0], h[0], lo[0] = 100.0, 100.5, 99.5

idx = pd.bdate_range("2020-01-01", periods=n)
ohlc = pd.DataFrame({"open": o, "high": h, "low": lo, "close": close,
                     "volume": 1e6}, index=idx).iloc[1:]

Planting a Known Premium

Now for the implied volatility. We set it to the volatility that will follow each day, plus exactly 4 points, plus a bit of noise so it isn't perfectly clean. If the pipeline pairs everything up correctly, the average premium it finds should come out close to +4.00.

IV_t = RV_{t \to t+h} + 0.04 + \varepsilon_t, \qquad \varepsilon_t \sim N(0,\, 0.015^2)

python

# IV = the realised vol that WILL follow, plus a 4-point premium, plus noise
TRUE_PREMIUM = 0.04
rv_fwd = forward_realised_vol(ohlc, HORIZON)
iv = (rv_fwd + TRUE_PREMIUM + rng.normal(0, 0.015, len(rv_fwd))).rename("iv")
iv = iv.ffill().bfill()

python

# IV = the realised vol that WILL follow, plus a 4-point premium, plus noise
TRUE_PREMIUM = 0.04
rv_fwd = forward_realised_vol(ohlc, HORIZON)
iv = (rv_fwd + TRUE_PREMIUM + rng.normal(0, 0.015, len(rv_fwd))).rename("iv")
iv = iv.ffill().bfill()

Checking It Comes Back Out

We run the fake data through exactly the same analyse() function the live version uses, then rebuild the premium series and check the average is within 1.2 points of the 4 we planted ( $|\overline{VRP} - 0.04| < 0.012$ ). The option being scored has an implied volatility of 40.45 against an at-the-money volatility of 38.00, so the skew adjustment gets exercised too.

python

print(report(analyse(iv, ohlc, market_iv=0.4045, atm_iv=0.3800)))

vrp_df = build_vrp_series(iv, ohlc, HORIZON)
recovered = vrp_baseline(vrp_df, 756, HORIZON)["mean"]
print(f"\n  injected premium   {100*TRUE_PREMIUM:+.2f} pts")
print(f"  recovered premium  {100*recovered:+.2f} pts")
assert abs(recovered - TRUE_PREMIUM) < 0.012, "pipeline failed to recover VRP"
print("  PASS — forward alignment is correct\n")

python

print(report(analyse(iv, ohlc, market_iv=0.4045, atm_iv=0.3800)))

vrp_df = build_vrp_series(iv, ohlc, HORIZON)
recovered = vrp_baseline(vrp_df, 756, HORIZON)["mean"]
print(f"\n  injected premium   {100*TRUE_PREMIUM:+.2f} pts")
print(f"  recovered premium  {100*recovered:+.2f} pts")
assert abs(recovered - TRUE_PREMIUM) < 0.012, "pipeline failed to recover VRP"
print("  PASS — forward alignment is correct\n")

The Earnings Round Trip

For the earnings maths we don't need a simulation at all - we can go backwards. We pick a true everyday volatility (35%) and a true earnings jump (6%), use the formula from earlier to work out what the implied volatilities of a 30-day and a 90-day option would be:

\sigma_f = \sqrt{0.35^2 + \frac{0.06^2}{30/365}}, \qquad \sigma_b = \sqrt{0.35^2 + \frac{0.06^2}{90/365}}

Then we hand those to the real functions to see if they get back 35% and 6%. Since there's no noise, they should match to 9 decimal places.

We also check the guard rails: a term structure with no earnings hump has to be refused rather than inventing a jump, and the expiry helpers are tested against a fixed date so the test doesn't start failing as time passes.

python

SIGMA_BASE, TRUE_J = 0.35, 0.06
T_f, T_b = 30 / 365, 90 / 365
iv_f = np.sqrt(SIGMA_BASE**2 + TRUE_J**2 / T_f)
iv_b = np.sqrt(SIGMA_BASE**2 + TRUE_J**2 / T_b)
rec_J = implied_jump_from_term_structure(iv_f, T_f, iv_b, T_b)
rec_base = strip_event_variance(iv_f, T_f, rec_J)
assert abs(rec_J - TRUE_J) < 1e-9, "jump identification failed"
assert abs(rec_base - SIGMA_BASE) < 1e-9, "event stripping failed"

# no event hump (back IV above front) must refuse, not invent a jump
try:
    implied_jump_from_term_structure(0.30, T_f, 0.35, T_b)
    raise AssertionError("should have refused a non-hump term structure")
except ValueError:
    pass

# expiry helpers, pinned to a fixed date so this stays deterministic
asof = date(2026, 9, 3)
exps = ["20260918", "20261016", "20270115", "20280121"]
assert abs(years_to_expiry("20270115", asof) - 134 / 365) < 1e-12
assert pick_back_expiry(exps, "20260918", asof=asof) == "20261016"
assert pick_back_expiry(exps, "20280121", asof=asof) is None

python

SIGMA_BASE, TRUE_J = 0.35, 0.06
T_f, T_b = 30 / 365, 90 / 365
iv_f = np.sqrt(SIGMA_BASE**2 + TRUE_J**2 / T_f)
iv_b = np.sqrt(SIGMA_BASE**2 + TRUE_J**2 / T_b)
rec_J = implied_jump_from_term_structure(iv_f, T_f, iv_b, T_b)
rec_base = strip_event_variance(iv_f, T_f, rec_J)
assert abs(rec_J - TRUE_J) < 1e-9, "jump identification failed"
assert abs(rec_base - SIGMA_BASE) < 1e-9, "event stripping failed"

# no event hump (back IV above front) must refuse, not invent a jump
try:
    implied_jump_from_term_structure(0.30, T_f, 0.35, T_b)
    raise AssertionError("should have refused a non-hump term structure")
except ValueError:
    pass

# expiry helpers, pinned to a fixed date so this stays deterministic
asof = date(2026, 9, 3)
exps = ["20260918", "20261016", "20270115", "20280121"]
assert abs(years_to_expiry("20270115", asof) - 134 / 365) < 1e-12
assert pick_back_expiry(exps, "20260918", asof=asof) == "20261016"
assert pick_back_expiry(exps, "20280121", asof=asof) is None

Running it prints the full decision row for the fake option, followed by the checks (I've trimmed a few printed notes from the output):

python

.venv/bin/python run_analysis.py --selftest

Output

==============================================================
  VARIANCE RISK PREMIUM — DECISION ROW
==============================================================
  Market IV           40.45 vol
  RV forecast         44.35 vol
  Gap                 -3.90 pts
  Typical gap          6.32 pts  (+/- 0.09 at 2 s.e.)
  Edge               -10.22 pts   z = -7.24
  ---> Underpaid to sell / relatively cheap to buy

  skew spread carried   +2.45 vol pts
  EWMA cross-check      34.81 vol

  Baseline (trailing 3y, ATM):
    observations        756
    mean premium        +3.87 pts
    median              +3.88 pts
    std dev             1.41 pts
    Newey-West s.e.     0.04 pts
    share positive      99.9%
    p10 / p90           +2.09 / +5.73
==============================================================

  injected premium   +4.00 pts
  recovered premium  +3.87 pts
  PASS — forward alignment is correct

  same data, BACKWARD-aligned: +4.00 pts

  implied earnings move  5.82%
  event-stripped vol     34.99%  (vs 40.45 raw)

  round trip: injected J 6.00% -> recovered 6.00%
              injected base 35.00 -> recovered 35.00
  PASS — decomposition inverts, guards hold, expiry helpers agree

python

.venv/bin/python run_analysis.py --selftest

Output

==============================================================
  VARIANCE RISK PREMIUM — DECISION ROW
==============================================================
  Market IV           40.45 vol
  RV forecast         44.35 vol
  Gap                 -3.90 pts
  Typical gap          6.32 pts  (+/- 0.09 at 2 s.e.)
  Edge               -10.22 pts   z = -7.24
  ---> Underpaid to sell / relatively cheap to buy

  skew spread carried   +2.45 vol pts
  EWMA cross-check      34.81 vol

  Baseline (trailing 3y, ATM):
    observations        756
    mean premium        +3.87 pts
    median              +3.88 pts
    std dev             1.41 pts
    Newey-West s.e.     0.04 pts
    share positive      99.9%
    p10 / p90           +2.09 / +5.73
==============================================================

  injected premium   +4.00 pts
  recovered premium  +3.87 pts
  PASS — forward alignment is correct

  same data, BACKWARD-aligned: +4.00 pts

  implied earnings move  5.82%
  event-stripped vol     34.99%  (vs 40.45 raw)

  round trip: injected J 6.00% -> recovered 6.00%
              injected base 35.00 -> recovered 35.00
  PASS — decomposition inverts, guards hold, expiry helpers agree

The earnings test is a genuine one - it starts from known true values and recovers them to 9 decimal places. The premium test is weaker than it looks, though. The synthetic implied volatility is built from the same forward-volatility function that the pipeline later subtracts, so the two cancel out: when I deliberately broke the forward alignment, shifted it by 5 days, or doubled the Yang-Zhang output, the test still passed. You can see a hint of this in the output - the backward-aligned version gives +4.00 too. It does catch sign errors, scaling errors and mistakes in the earnings maths, but a better test would generate the implied volatility from the simulation's true volatility instead.

The Newey-West error bar here (0.04 points) is also unrealistically tight for the same reason - on real data the band is far wider.

Running It on a Real Option

To see what happens on real data, I ran it on an Apple put: the $340 strike expiring on 16 October 2026, which was the monthly expiry closest to a month away. Apple had closed at $341.07, so this is basically at the money. I ran it on 27 September 2026 using IBKR's delayed data, since I don't have a live options data subscription - so the implied volatility is a delayed quote rather than a live one, and the tool says so.

python

.venv/bin/python run_analysis.py AAPL 20261016 340 P

Output

NOTE: market data type 4 — this IV may be delayed or a stale close, not a live quote. Do not trade off it without confirming.

==============================================================
  VARIANCE RISK PREMIUM — DECISION ROW
==============================================================
  Market IV           20.84 vol
  RV forecast         24.84 vol
  Gap                 -4.00 pts
  Typical gap         -2.06 pts  (+/- 2.45 at 2 s.e.)
  Edge                -1.94 pts   z = -0.22
  ---> NO EDGE — inside the noise band of the baseline itself

  skew spread carried   +0.00 vol pts
  EWMA cross-check      22.76 vol

  Baseline (trailing 3y, ATM):
    observations        756
    mean premium        -2.06 pts
    median              +0.05 pts
    std dev             8.90 pts
    Newey-West s.e.     1.22 pts
    share positive      50.3%
    p10 / p90           -9.39 / +5.21

  Earnings decomposition: NOT APPLIED
    term structure not in event-hump configuration
    IV below is raw and may contain an event jump.
==============================================================

python

.venv/bin/python run_analysis.py AAPL 20261016 340 P

Output

NOTE: market data type 4 — this IV may be delayed or a stale close, not a live quote. Do not trade off it without confirming.

==============================================================
  VARIANCE RISK PREMIUM — DECISION ROW
==============================================================
  Market IV           20.84 vol
  RV forecast         24.84 vol
  Gap                 -4.00 pts
  Typical gap         -2.06 pts  (+/- 2.45 at 2 s.e.)
  Edge                -1.94 pts   z = -0.22
  ---> NO EDGE — inside the noise band of the baseline itself

  skew spread carried   +0.00 vol pts
  EWMA cross-check      22.76 vol

  Baseline (trailing 3y, ATM):
    observations        756
    mean premium        -2.06 pts
    median              +0.05 pts
    std dev             8.90 pts
    Newey-West s.e.     1.22 pts
    share positive      50.3%
    p10 / p90           -9.39 / +5.21

  Earnings decomposition: NOT APPLIED
    term structure not in event-hump configuration
    IV below is raw and may contain an event jump.
==============================================================

The verdict is "no edge": the option is 1.94 points cheaper than usual, but that's well inside the error bar of ±2.45 (Newey-West standard error rounded down to 1.22, and twice that gives 2.45 rounding up). Notice how much wider that band is than the 0.04 in the self-test - real data is a lot messier.

The most interesting part is the baseline. You'd expect the premium to be positive, since options are insurance, but Apple's average came out negative (-2.06 points), even though the median is basically zero (+0.05).

Sorting the premium series shows why: the worst readings are all from late March 2025, when Apple's implied volatility was about 23% but the stock then moved at about 75% volatility over the following month (the April 2025 tariff crash). That's around -50 points for a few weeks straight, which is enough to drag a 3-year average below zero on its own. This is exactly what selling insurance looks like - you collect a small premium most of the time, and then occasionally you pay out for a disaster. It's also why the median and the p10/p90 range are worth looking at, not just the mean.

There are a couple of other things to notice. The skew spread is zero because $340 is the at-the-money strike itself, so there's no strike adjustment to make. And the earnings decomposition wasn't applied: the October option's implied volatility wasn't higher than the longer-dated (6 November) one's, so the market isn't pricing an earnings jump into this contract.

A Stock That Missed the Crash: McDonald's

Since one month of 2025 was doing so much damage to Apple's baseline, I wanted to see what happens on a stock that wasn't hit as hard. Defensive stocks (things people keep buying in a downturn, like food, household goods and electricity) tend to fall less in a crash, so I ran it on an at-the-money McDonald's put: the $237.50 strike, with the same 16 October 2026 expiry. McDonald's had closed at $236.50.

python

.venv/bin/python run_analysis.py MCD 20261016 237.5 P

Output

NOTE: market data type 4 — this IV may be delayed or a stale close, not a live quote. Do not trade off it without confirming.

==============================================================
  VARIANCE RISK PREMIUM — DECISION ROW
==============================================================
  Market IV           19.70 vol
  RV forecast         18.74 vol
  Gap                  0.96 pts
  Typical gap         -0.33 pts  (+/- 1.03 at 2 s.e.)
  Edge                 1.29 pts   z = +0.33
  ---> Overpaid to sell vol relative to the historical premium

  skew spread carried   +0.00 vol pts
  EWMA cross-check      24.75 vol

  Baseline (trailing 3y, ATM):
    observations        756
    mean premium        -0.33 pts
    median              +0.17 pts
    std dev             3.89 pts
    Newey-West s.e.     0.51 pts
    share positive      52.6%
    p10 / p90           -4.82 / +3.82

  Earnings decomposition: NOT APPLIED
    no longer-dated expiry available to identify J.
    IV below is raw and may contain an event jump.
==============================================================

python

.venv/bin/python run_analysis.py MCD 20261016 237.5 P

Output

NOTE: market data type 4 — this IV may be delayed or a stale close, not a live quote. Do not trade off it without confirming.

==============================================================
  VARIANCE RISK PREMIUM — DECISION ROW
==============================================================
  Market IV           19.70 vol
  RV forecast         18.74 vol
  Gap                  0.96 pts
  Typical gap         -0.33 pts  (+/- 1.03 at 2 s.e.)
  Edge                 1.29 pts   z = +0.33
  ---> Overpaid to sell vol relative to the historical premium

  skew spread carried   +0.00 vol pts
  EWMA cross-check      24.75 vol

  Baseline (trailing 3y, ATM):
    observations        756
    mean premium        -0.33 pts
    median              +0.17 pts
    std dev             3.89 pts
    Newey-West s.e.     0.51 pts
    share positive      52.6%
    p10 / p90           -4.82 / +3.82

  Earnings decomposition: NOT APPLIED
    no longer-dated expiry available to identify J.
    IV below is raw and may contain an event jump.
==============================================================

This time the verdict is "overpaid to sell": the option is 1.29 points more expensive than usual, which just clears the error bar of ±1.03. The band is much tighter than Apple's because McDonald's premium barely moves around - its standard deviation is 3.89 points compared to Apple's 8.90. It's still a small effect though (z = +0.33), so I wouldn't read too much into it.

The earnings decomposition wasn't applied here either, but for a different reason - the tool couldn't find a longer-dated expiry to compare against, which is strange for a stock with plenty of listed options. McDonald's usually reports earnings after mid-October, so it probably doesn't affect this result, but it's something I need to look into.

Before picking McDonald's, I checked how much March and April 2025 hurt the premium for a few defensive stocks, with Apple for comparison (all in volatility points, over the same 3 years the tool uses):

Stock	3-year mean	Median	Worst reading (date)	Mar–Apr 2025 average
McDonald’s (MCD)	−0.33	+0.17	−13.70 (27 Sep 2024)	−2.03
Duke Energy (DUK)	−0.45	−0.12	−13.38 (24 Mar 2025)	−4.08
Procter & Gamble (PG)	−0.41	−0.11	−13.30 (26 Mar 2025)	−4.59
Coca-Cola (KO)	−0.99	−0.29	−17.43 (26 Mar 2025)	−5.88
Verizon (VZ)	−1.74	−0.64	−21.43 (14 Feb 2025)	−6.84
Johnson & Johnson (JNJ)	−0.41	−0.12	−22.77 (24 Mar 2025)	−9.39
Apple (AAPL)	−2.06	+0.05	−51.90 (25 Mar 2025)	−20.86

All of them were hit much less than Apple, and McDonald's was the least affected - its worst reading isn't even from 2025. However, something else stands out: even for these stocks, the average premium is around zero, and only about half of the readings are positive. Since selling options should normally earn a clearly positive premium, the 2025 crash can't be the whole story. My guess is that something systematic is going on, either in IBKR's implied volatility series or in how it compares to the Yang-Zhang estimate, and rebuilding the baseline with a different measure of realised volatility would be a good way to find out.

Limitations

There are a few things this tool can't tell you:

It checks the level, not the shape. Today's extra premium for being away from the money is taken at face value and never compared against its own history. If the skew is unusually steep right now, the tool treats it as normal. So the answers are sharpest near the money and weakest far away from it. Fixing this needs historical data for every individual strike, which IBKR doesn't keep - it's a data purchase rather than a code change.
The earnings clean-up only happens on one side. Today's price has the earnings bump removed, but the historical premium still has earnings periods mixed into it, since stripping those would need a history of past earnings dates. It's right in direction, but not exact.
The horizon is always 21 trading days. The forecast and the premium are both one-month numbers, whatever the option's expiry, so a long-dated option is being compared against the wrong horizon. This barely matters for options expiring in about a month, but it matters a lot for ones expiring a year out.
The baseline premium is suspiciously low. Selling options should earn a clearly positive premium on average, but every stock I checked came out around zero, and Apple's was negative. The 2025 crash explains part of this, but not all of it, so something systematic may be off - either in IBKR's implied volatility series or in how it lines up with the Yang-Zhang estimate. Until that's understood, the typical gap (and so the edge) could be shifted for every stock.
No backtest. IBKR doesn't provide historical data for expired option contracts, so there's no way to check how the verdicts would have performed in the past with this data source.

One practical note: the live output also reports the gap between the bid and ask implied volatilities. If that's more than about 5 points, there isn't really a meaningful price to analyse - the tool warns you about this, but doesn't stop you.

Thanks for reading - I hope this helps, and don't forget to check out the full source code on GitHub too!