In this project we will get more into volatility and try to use statistical models like GARCH and HAR-RV to forecast realized volatility and then work with strategies and signals that come from the difference between market implied volatility and realized volatility.
Introduction to GARCH and implementation
rolling forecasts using GARCH, GJR GARCH and EGARCH models
Implamantation of HAR-RV
Computation of Variance risk premium and interpretation
using VRP as a trading signal for straddle strategies and volatility trading
In Project 4, the central object was the option price: given an option chain, rates, forwards, maturities, and a pricing model, we solved for implied volatility and Greeks. Here the direction becomes more trading-oriented: we ask whether the volatility implied by option prices is high or low relative to a model-based forecast of future realized variance, and then we translate that signal into simple straddle overlays.
This is why the models are evaluated mainly in variance space. It is also why the VRP is defined as an implied-variance spread rather than only as a volatility spread. Volatility is easier to visualize, but variance is the natural object for multi-period forecasting and option-risk-premium logic.
physical measure
There is also a separation between the physical measure and the option-implied measure. GARCH and HAR-RV are estimated from realized return histories, so their forecasts are physical forecasts under the historical data process. Implied volatility comes from market prices and includes risk premia, supply-demand pressure, jump fear, margin effects, and hedging flows. Therefore, IV can be high or low relative to a disciplined physical forecast, and this comes from the difference between theory and what the volatility can be in a stationary world and what it is in real world markets.
1) Data
We use the same option chain, treasury risk free rate and historical prices of SPX in this project so we can have access to implied volatility solved in project 4. We don’t need SPY because we are not going to deal with greeks and hedging with an underlying in this notebook. you can donload the options chain from here: SPX EOD option chains from 2022 and 2023 and risk free rate from here: treasury par yields and SPX historical data is from yfinance API
1.1 building future realized volatility
for evaluating the forecasts in the next sections, we have to construct the targets and the real volatility of each time. since we want to have forecasts for different horizons, we have to find the realized volatility of all the horizons for each date in our 2 year window.
The key transformation is from price levels to log returns:
\[
r_t=\log S_t-\log S_{t-1}.
\]
For each horizon \(h\in\{1,5,10,21,42,63\}\), the code builds the future realized variance target
Later, we clarify why we use squared returns as the estimator of future variance.
This is a forward target, not a realized volatility estimate. That means date \(t\) receives the realized variance that will be observed over the next \(h\) trading days. this is the value that the forecasting model is trying to predict.
The trailing measure is observable at date \(t\) and can be used as a feature. The forward measure is not observable at date \(t\) and is only used for validation. Plotting them together gives a visual sense of how hard the forecasting problem is: trailing volatility reacts to information already seen, while forward volatility can jump if the next month contains a shock.
Now we repeat the process of solving the implied volatility from Black Scholes formula and we use our library quantfinlab and use the Let’s Be Rational inspired solver that we built in notebook 4.
\[
C_{\text{model}}(\sigma) = C_{\text{mkt}}
\]
In the quotes cleaning process, We focus on near-ATM IV because we want a clean volatility comparison, not a full smile model. Deep OTM options include strong jump risk and crash insurance premia. Deep ITM options often have worse liquidity and wider microstructure noise. Near-ATM options are usually the most direct way to extract the market’s central volatility expectation.
A volatility forecast can be very sophisticated, but if the implied-volatility input is noisy, stale, or taken from illiquid options, the signals made from it become unreliable.
We now construct and fit a GARCH model manually before using the arch package. This is not because the manual implementation is faster or more complete. It is because showing a transparent and clear implementation of model instead of just using a black box function from this package. After this we use arch package for all the other models and forecasts.
Daily asset returns exhibit:
Zero autocorrelation in levels (weak efficiency, today’s return is not correlated to tommorow or yesterday).
Volatility clustering: large (small) returns tend to be followed by large (small) returns of either sign.
Heavy tails distribution (excess kurtosis).
Slow decay of autocorrelation of squared returns (long memory in volatility).
3.1 Conditional variance
If \(\{r_t\}_{t \in \mathbb{Z}}\) is a return process adapted to the filtration \(\mathcal{F}_{t-1}\) (all information up to time \(t-1\)).
For a random variable \(r_t\) and information set \(\mathcal{F}_{t-1}\):
This is why we used squared returns as the estimator for future realized volatility and as the target of GARCH models. Under the assumption that return has mean of 0 and is unpredictable.
Then we can get to The unconditional variance: \(\operatorname{Var}(r_t) = \mathbb{E}[r_t^2]\).
And then We can define the conditional variance as:
where \(\{z_t\}\) is an i.i.d. sequence with zero mean and unit variance (often taken as standard normal or Student‑\(t\)). It is a random shock added to volatility that we suppose will construct the unpredictable return.
\(\sigma_t\) is \(\mathcal{F}_{t-1}\)-measurable (predictable). Then:
But \(\mathbb{E}[r_{t-1}^2] = \mathbb{E}[\sigma_{t-1}^2]\) (since \(z_{t-1}\) has unit variance and is independent). Let \(\bar{\sigma}^2 = \mathbb{E}[\sigma_t^2]\) (constant under stationarity). Then:
which is positive only if \(1 - \alpha - \beta > 0\) or the rearranged constraint \(\alpha + \beta < 1\).
3.4 The ARMA representation for squared returns
If \(\nu_t = r_t^2 - \sigma_t^2\). Then \(\nu_t\) is a zero‑mean white noise (though not i.i.d. – it has heteroskedasticity). Substituting \(\sigma_t^2 = r_t^2 - \nu_t\) into the GARCH(1,1) equation:
This is an ARMA(1,1) model for \(r_t^2\) with AR (Autoregressive) coefficient \(\phi = \alpha+\beta\) and MA (Moving average) coefficient \(\theta = -\beta\).
The persistence of squared returns is governed by \(\alpha+\beta\). if it is close to 1, shocks to \(r_t^2\) decay very slowly, matching the stylised fact of long memory in volatility.
3.5 Maximum likelihood estimation (Gaussian)
Given observed returns \(r_1, \dots, r_T\), we estimate the parameters set \(\theta = (\omega, \alpha, \beta)\) by maximising the log‑likelihood.
For \(z_t \stackrel{\text{i.i.d.}}{\sim} \mathcal{N}(0,1)\):
To compute \(\sigma_t^2\) for \(t=1\), we need \(\sigma_0^2\) (and \(r_0^2\), but that can be set to its sample mean).
A common method is backcasting: We can choose \(\sigma_0^2\) as the exponentially weighted moving average of early squared returns:
where \(\lambda\) is a decay factor (usually 0.94 for daily data) and \(m\) is a moderate number (like 100). This gives more weight to recent observations before the sample start.
3.7 Multi‑step variance forecasting
We need forecasts of cumulative variance over horizons \(h = 1,5,10,21,42,63\) days.
The one step ahead forecast is directly computed from:
Annualized volatility would be \(\sqrt{252 \cdot \widehat{RV}_{t,1}}\) for one day ahead, but for multi‑day horizons we use the sum directly.
Show code
returns_series = df_market.loc[df_market["spx_ret"].notna(), ["trade_date", "spx_ret"]].copy().reset_index(drop=True)ret_pct =100.0* returns_series["spx_ret"].to_numpy(dtype=float)def garch_backcast(resids, decay=0.94, horizon=75): x = np.asarray(resids, dtype=float) **2 n =min(len(x), horizon)if n ==0:return1.0 w = decay ** np.arange(n) w = w / w.sum()returnfloat(np.sum(w * x[:n][::-1]))def garch11_filter_backcast(returns_pct, omega, alpha, beta, backcast=None): eps = np.asarray(returns_pct, dtype=float) n =len(eps) sigma2 = np.empty(n, dtype=float)if backcast isNone: backcast = garch_backcast(eps) sigma2[0] =max(float(backcast), 1e-8)for t inrange(1, n): sigma2[t] = omega + alpha * eps[t -1] **2+ beta * sigma2[t -1] sigma2[t] =max(sigma2[t], 1e-10)return sigma2def garch11_negloglik_normal(params, returns_pct, backcast=None): omega, alpha, beta = [float(x) for x in params]if omega <=0.0or alpha <0.0or beta <0.0or alpha + beta >=0.999:return1e16 sigma2 = garch11_filter_backcast(returns_pct, omega, alpha, beta, backcast=backcast) ll =-0.5* (np.log(2.0* np.pi) + np.log(sigma2) + (returns_pct **2) / sigma2) value =-float(np.sum(ll))return value if np.isfinite(value) else1e16def garch11_multi_horizon_variance(last_sigma2, last_eps, omega, alpha, beta, horizons): horizons =sorted(set(int(h) for h in horizons)) out = {} h1 = omega + alpha *float(last_eps) **2+ beta *float(last_sigma2) out[1] = h1 prev = h1for h inrange(2, max(horizons) +1): prev = omega + (alpha + beta) * prev out[h] = prevreturn {h: out[h] for h in horizons}def fit_manual_garch11_normal(returns_pct): y = np.asarray(returns_pct, dtype=float) y = y[np.isfinite(y)] backcast = garch_backcast(y) var0 =max(float(np.var(y, ddof=1)), 1e-6) starts = []for alpha0, beta0 in [(0.05, 0.90), (0.08, 0.90), (0.10, 0.85), (0.03, 0.95), (0.12, 0.80)]:if alpha0 + beta0 <0.999: starts.append(np.array([max((1.0- alpha0 - beta0) * var0, 1e-8), alpha0, beta0], dtype=float)) bounds = [(1e-12, 100.0* var0), (1e-8, 0.999), (1e-8, 0.999)] cons = [{"type": "ineq", "fun": lambda x: 0.999- x[1] - x[2]}] best =Nonefor x0 in starts:try: res = minimize( garch11_negloglik_normal, x0=x0, args=(y, backcast), method="SLSQP", bounds=bounds, constraints=cons, options={"maxiter": 2000, "ftol": 1e-12})if best isNoneor res.fun < best.fun: best = resexceptException:passif best isNone:raiseRuntimeError("manual GARCH fit failed") omega, alpha, beta = [float(v) for v in best.x] sigma2 = garch11_filter_backcast(y, omega, alpha, beta, backcast=backcast) std_resid = y / np.sqrt(np.clip(sigma2, 1e-12, None)) ll =-garch11_negloglik_normal([omega, alpha, beta], y, backcast) n =len(y) k =3return {"success": bool(best.success), "message": str(best.message),"omega": omega, "alpha": alpha, "beta": beta,"persistence": alpha + beta, "backcast": backcast,"sigma2": sigma2, "std_resid": std_resid, "loglik": ll,"aic": 2.0* k -2.0* ll, "bic": math.log(n) * k -2.0* ll,"unc_var": omega /max(1.0- alpha - beta, 1e-8)}manual_fit = fit_manual_garch11_normal(ret_pct)am_manual_check = arch_model(ret_pct, mean="Zero", vol="GARCH", p=1, q=1, dist="normal", rescale=False)res_manual_check = am_manual_check.fit(disp="off", show_warning=False)display(pd.DataFrame([ {"model": "manual_garch11_normal","success": manual_fit["success"],"omega": manual_fit["omega"],"alpha": manual_fit["alpha"],"beta": manual_fit["beta"],"persistence": manual_fit["persistence"],"unc_var_pct2": manual_fit["unc_var"],"aic": manual_fit["aic"],"bic": manual_fit["bic"]}, {"model": "arch_garch11_normal","success": True,"omega": float(res_manual_check.params.get("omega", np.nan)),"alpha": float(res_manual_check.params.get("alpha[1]", np.nan)),"beta": float(res_manual_check.params.get("beta[1]", np.nan)),"persistence": float(res_manual_check.params.get("alpha[1]", np.nan) + res_manual_check.params.get("beta[1]", np.nan)),"unc_var_pct2": float(res_manual_check.params.get("omega", np.nan) /max(1.0- res_manual_check.params.get("alpha[1]", np.nan) - res_manual_check.params.get("beta[1]", np.nan), 1e-8)),"aic": float(res_manual_check.aic),"bic": float(res_manual_check.bic)}]))
model
success
omega
alpha
beta
persistence
unc_var_pct2
aic
bic
0
manual_garch11_normal
True
0.010925
0.095796
0.899325
0.995121
2.239160
65,498.576815
65,522.914655
1
arch_garch11_normal
True
0.010923
0.095777
0.899349
0.995127
2.241470
65,498.126263
65,522.464103
We tried to minimize the negative log likelihood using a numerical optimizer (L‑BFGS‑B with bounds \([\epsilon, 1-\epsilon]\) for \(\alpha,\beta\)). As we can see our fitted model is pretty close to the implementation of arch package and got near exact results.
The fitted parameters are:
Parameter
Manual
arch
\(\omega\)
0.010925
0.010923
\(\alpha\)
0.095796
0.095777
\(\beta\)
0.899325
0.899349
\(\alpha\) is small (≈0.096): A large shock today increases tomorrow’s variance by \(\alpha\) times the squared shock.
\(\beta\) is large (≈0.899): Most of tomorrow’s variance is inherited from today’s variance (volatility is persistent).
The sum \(\alpha+\beta \approx 0.995\) which is less than 1 and implies that shocks take a long time to decay.
We can interpret it like this: A large negative or positive return creates a jump in \(r_t^2\), which enters through \(\alpha r_t^2\), and then the shock decays through the \(\beta\) channel over many days.
The absolute parameter differences are extremely small, and the correlation between our and arch implementation is 0.999991. The only visible difference is an initialization effect. After that, the two implementations are practically the same.
The residual histogram being close to normal in the center is expected, but the QQ plot reveals that the model might perform close with normality but it misses the tails. This is a common trap. A model can look visually good for most of the observations and still be dangerous for option trading if it doesn’t represent extreme moves as good. Options are nonlinear instruments, so tail errors matter more than they would in a simple linear regression.
The standardized residuals:
\[
z_t=\frac{r_t}{\widehat{\sigma}_t}
\]
We compare empirical residual quantiles with normal quantiles \[
Q_{\widehat{z}}(p)
\quad \text{versus} \quad
\Phi^{-1}(p).
\]
If the residuals were exactly standard normal, the points would lie on the line completely. The deviation in the tails means the standardized shocks still have excess kurtosis and asymmetry.This result can be a sign that we can use different distrinutions that can represent tails better than normal distribution. One of the main ones used for tail tracking is Student-t distribution
3.8 Beyond normality
A Student-\(t\) distribution with degrees of freedom \(\nu\) has heavier tails than the normal distribution. When \(\nu\) is finite and not very large, the model assigns more probability to large standardized shocks:
\[
P(|Z_t|>x)
\text{ is larger under Student-}t
\text{ than under } \mathcal{N}(0,1).
\]
This is one of the directions for equity-index return modeling that can perform better at tails than normal GARCH.
GARCH can explain volatility clustering, but it does not automatically explain tail shape. A return can be standardized by a strong volatility filter and still have extreme residual in the tails. That is the reason we can use Student-\(t\) here:
\[
z_t \sim t_{\nu},
\]
where lower \(\nu\) means heavier tails. For an equity index’s returns, we except much heavier tailed behavior than normal, but not so extreme that the variance becomes undefined.
While GARCH(1,1) captures most financial volatility dynamics, a GARCH(2,2) allows for two lags of both squared returns and past variances, might be better fitting assets with more complex memory structures:
Here \(\gamma>0\) means negative returns increase future variance more than positive returns of the same magnitude. This means the negative shocks that make the price of an asset fall have higher effect on variance than positive shocks that make the price rise. This is the standard leverage effect pattern in equity indices.
Because the model is written in \(\log\sigma_t^2\), we ensure the variance is automatically positive without requiring the same parameter restrictions as standard GARCH.
We also compare normal and Student-\(t\) distributions for all three of these models:
GARCH(1,1) Normal
GARCH(1,1) t
GARCH(2,2) Normal
GARCH(2,2) t
GJR-GARCH(1,1,1) Normal
GJR-GARCH(1,1,1) t
EGARCH(1,1) Normal
EGARCH(1,1) t
There are so many more GARCH variants that can be even more powerful and perform better in modeling, like Mixture GARCH models, multivariate, Markov-switching, CC_GARCH, DCC-GARCH etc, which have their own features and abilities. for this project we will focus on the main classic ones. We also just implemented GARCH(1,1) manualy for clarity of model and for the rest of the models, we use arch package.
GJR-GARCH(1,1) with Student-\(t\) has the best AIC and BIC.
The Student-\(t\) models strongly improve likelihood relative to the normal models.
The estimated GJR asymmetry, with \(\gamma_1\approx0.113\) in the best model.
The plain normal GARCH(1,1) is useful as a baseline, but it is not the best full-sample description.
The Student-\(t\) models estimate degree of freedom \(\nu\) around 5.8–6.2, which is strong evidence that the return distribution is much heavier tailed than Gaussian (normal).
and the ARCH-LM test checks whether squared residuals still forecast themselves. The p-values are mostly very small, especially for raw standardized residual autocorrelation. This is not surprising in a very long sample. With more than 20,000 daily observations, even small misspecifications become statistically significant.
the GJR Student-\(t\) model is the best full-sample likelihood model
asymmetry and fat tails matter
no single parametric volatility model removes all dependence completely
full sample fit is not the same as out of sample forecast quality. because of this, we will get to rolling forecasts in the next sections to be able to compare models in real condition.
Persistent is different for each model type and have to be interpreted different. The stronger interpretation is that asymmetric and fat-tailed models describe the data better than the plain normal GARCH baseline.
As we can see from the plots, different models often agree on the dates of volatility spikes but differ in the height and decay of those spikes. For ana pplication like option trading, this matters because a small difference in annualized volatility becomes a larger difference in variance:
The best full sample model is GJR-GARCH with Student-\(t\) distribution. That means two empirical features matter at the daily SPX level:
asymmetric response to negative shocks + fat-tailed standardized shocks
5) HAR-RV model
We now add a non-GARCH model. HAR-RV is the heterogeneous autoregressive model of realized volatility. GARCH models conditional variance as a long run mean reverting process driven by shocks. HAR-RV models realized variance directly using multiple realized volatility horizons.
The whole HAR idea is that volatility has heterogeneous memory. Short horizon traders react to daily volatility, medium horizon participants react to weekly volatility, and slower allocators react to monthly volatility. A simple HAR-RV regression uses
There is also an important detail for preventing leaking future data to model. For a horizon-\(h\) forecast made at date \(t\), the target \(RV_{s,s+h}^{\text{sum}}\) for a training date \(s\) is only known after \(s+h\). Therefore, when forecasting at \(t\), the training target cannot include dates whose realized window contains the unknown future, so the regression only trains on targets that would already have been fully realized. That detail is easy and obvious, but it is essential. Without it, the HAR model would accidentally use future realized variance information and would look better than it should in a real time trading setting.
We use HAR-RV along GARCH models because it’s intentionally simple. It doesn’t estimate a nonlinear variance recursion and it doesn’t require a distributional assumption for returns. It says that future realized variance can be forecast from realized variance measured at different past horizons.
HAR-RV is also naturally connected to market behavior. Daily volatility may capture very recent shocks. Weekly volatility may capture short risk management cycles. Monthly volatility may capture slower repricing and institutional allocation changes. The model is linear, but its features encode a multi scale view of volatility.
HAR-type models are popular in realized-volatility forecasting because they are simple, interpretable, and often surprisingly competitive.
Show code
def make_har_features_visible(rv_daily, weekly_window=5, monthly_window=22, use_log=True, eps=eps): rv = pd.Series(rv_daily).copy() rv.index = pd.to_datetime(rv.index, errors="coerce") rv = pd.to_numeric(rv, errors="coerce").replace([np.inf, -np.inf], np.nan).sort_index() features = pd.DataFrame(index=rv.index) features["rv_daily"] = rv features["rv_weekly"] = rv.rolling(weekly_window, min_periods=weekly_window).mean() features["rv_monthly"] = rv.rolling(monthly_window, min_periods=monthly_window).mean()if use_log: features = np.log(features.clip(lower=eps)) features = features.rename(columns=lambda c: "log_"+ c)return featuresdef forward_realized_var_sum(returns, horizon): r = pd.Series(returns).copy() vals = pd.to_numeric(r, errors="coerce").to_numpy(dtype=float) out = np.full(len(vals), np.nan, dtype=float) h =int(horizon)for i inrange(len(vals) - h): window = vals[i +1 : i +1+ h]if np.isfinite(window).all(): out[i] =float(np.sum(window * window))return pd.Series(out, index=r.index, name=f"realized_var_sum_{h}")def fit_har_ols_visible(features, target_sum, use_log=True, eps=eps): data = features.join(target_sum.rename("target"), how="inner").replace([np.inf, -np.inf], np.nan).dropna()iflen(data) < features.shape[1] +20:raiseValueError("not enough HAR-RV observations") y = pd.to_numeric(data["target"], errors="coerce")if use_log: y = np.log(y.clip(lower=eps)) X = np.column_stack([np.ones(len(data)), data[features.columns].to_numpy(dtype=float)]) beta, *_ = np.linalg.lstsq(X, np.asarray(y, dtype=float), rcond=None)return pd.Series(beta, index=["const", *features.columns])def predict_har_sum_visible(params, feature_row, use_log=True, eps=eps): x = np.r_[1.0, feature_row.to_numpy(dtype=float)] pred =float(np.dot(params.to_numpy(dtype=float), x))if use_log: pred =float(np.exp(pred))returnmax(pred, eps)def rolling_har_forecasts_visible(returns_df, signal_dates, horizons=horizons, train_window=756, annualization=annualization, use_log=True, eps=eps): ser = returns_df.set_index("trade_date")["spx_ret"].dropna().sort_index() rv_daily = ser * ser features = make_har_features_visible(rv_daily, use_log=use_log, eps=eps) targets = {int(h): forward_realized_var_sum(ser, int(h)) for h in horizons} dates = pd.DatetimeIndex(ser.index) date_to_pos = {pd.Timestamp(d).normalize(): i for i, d inenumerate(dates)} records = []for signal_date in pd.DatetimeIndex(pd.to_datetime(signal_dates)).normalize().unique().sort_values(): pos = date_to_pos.get(pd.Timestamp(signal_date).normalize())if pos isNoneor pos < train_window or pos +max(horizons) >=len(ser):continueif signal_date notin features.index or features.loc[[signal_date]].isna().any(axis=None):continuefor h in horizons: h =int(h) known_end_pos = pos - hif known_end_pos <=0:continue train_start =max(0, known_end_pos -int(train_window) +1) train_dates = dates[train_start : known_end_pos +1]try: params = fit_har_ols_visible(features.loc[train_dates], targets[h].loc[train_dates], use_log=use_log, eps=eps) forecast_sum = predict_har_sum_visible(params, features.loc[signal_date], use_log=use_log, eps=eps)exceptException:continue realized_sum =float(targets[h].loc[signal_date]) forecast_ann = annualization / h * forecast_sum realized_ann = annualization / h * realized_sum if np.isfinite(realized_sum) else np.nan records.append({"date": signal_date,"model": "har_rv","horizon": h,"forecast_var_sum": forecast_sum,"forecast_var_ann": forecast_ann,"forecast_vol_ann": float(np.sqrt(max(forecast_ann, 0.0))),"realized_var_sum": realized_sum,"realized_var_ann": realized_ann,"realized_vol_ann": float(np.sqrt(max(realized_ann, 0.0))) if np.isfinite(realized_ann) else np.nan,"n_train": int(len(train_dates))})return pd.DataFrame(records).sort_values(["date", "horizon"]).reset_index(drop=True)har_feature_preview = make_har_features_visible(df_market.set_index("trade_date")["spx_ret"].dropna() **2).dropna().head()
6) rolling forecasts
For every signal date and every forecast horizon, we refit models using a rolling training window of 756 trading days, approximately three years:
\[
\mathcal{T}_t \{t-755,\ldots,t\}.
\]
The GARCH models produce a path of one step ahead conditional variances:
Rolling refits are expensive but necessary. thousands of GARCH and HAR-RV estimations are computationally heavy because we refit many models many times. A full-sample fit would be faster, but it would not represent the information set of a trader on each option date.
The structure is:
date t –> fit on past 756 returns –> forecast horizons –> store realized target for later evaluation
This gives a realistic forecast panel but introduces two practical complications.
A full-sample model can be useful for education, but it is not a trading test. In a trading test, the model parameters used on date \(t\) must be estimated using only information available before date \(t\). The rolling window creates exactly that discipline:
\[
\{\text{returns up to }t\}
\longrightarrow
\widehat{\theta}_t
\longrightarrow
\widehat{RV}_{t,t+h}.
\]
The 756-day training window is roughly three trading years. That is long enough to estimate GARCH-family parameters with some stability, but short enough to adapt when volatility regimes change.
The signal dates are the dates for which the notebook has a clean near-30-day ATM IV observation. In this sample, there are 505 signal dates from 2022-01-03 to 2023-12-29, which is the wondow that we have option chain data for.
On 2022-01-03, very short realized variance over the next day is tiny, but longer-horizon realized volatility is meaningfully higher. This shows why single horizon forecast can be misleading. A model may look wrong at \(h=1\) because the next day is quiet, while still being directionally reasonable over a month.
We have forecasted volatility using 9 models. Eight GARCH models and one HAR-RV. We can see a number of 26,784 forecast rows. This is slightly less than \(505\times9\times6 = 27,270\) because the end of the sample cannot support all forward realized windows up to 63 days. those final dates would not have a complete future target.
7) Forecast evaluation: QLIKE, Mincer–Zarnowitz calibration, and Diebold–Mariano
Now we want to evaluate the models under each horizon and find the best model for each one.
where \(\widehat{RV}_t\) is the forecasted variance and \(RV_t\) is the target realized volatility.
Lower QLIKE is better. QLIKE is very useful for variance forecasts because it penalizes under forecasting large realized variance very strongly while remaining appropriate when realized variance is a noisy representive of variance. We compute QLIKE on summed variance over the horizon, not directly on annualized volatility.
A perfectly calibrated unbiased forecast would have
\[
\alpha=0,\qquad \beta=1.
\]
Which means forecast and realized variance are exactly the same. In practice, \(\beta<1\) often means the forecast is too variable or too aggressive, while \(\beta>1\) often means the forecast underreacts to changes in realized variance.
7.3 Diebold–Mariano
We use this method directly for comparing models together. The Diebold–Mariano comparison is based on a loss differential:
\[
d_t=L_{t}^{(a)}-L_{t}^{(b)}.
\]
If the mean of \(d_t\) is significantly different from zero, the two forecasts have statistically different predictive loss. We use a long-run variance estimate because multi horizon forecast losses are serially dependent.
First, by average rank across the displayed metrics and horizons, EGARCH normal has the best overall rank, followed by HAR-RV. This does not mean EGARCH normal wins every horizon. It means it is best overall if we combine all the scores.
Second, at the 21-day horizon, GJR normal has the best QLIKE, with QLIKE about -4.933, while EGARCH normal has the lowest RMSE variance error among the top group. We can also see meaningful calibration differences from 21-day table. For example, EGARCH normal has a Mincer–Zarnowitz slope around 0.819, while HAR-RV has a slope above 1.20. This suggests that HAR-RV is less aggressive in some regimes and needs a larger slope to explain realized variance.
We can see the same story from heatmap. Shorter horizons are dominated by asymmetric GARCH-style models, while longer horizons become more mixed and HAR-RV becomes competitive. GARCH reacts quickly to return shocks, while HAR-RV often performs well at medium horizons because it combines daily, weekly, and monthly memory.
The QLIKE values in the tables are negative, which may look odd if you expect a loss to be positive. There is no problem here. QLIKE contains \(\log(\widehat{RV})\), and realized variances are often much smaller than 1 when returns are in decimal units. Therefore \(\log(\widehat{RV})\) can be strongly negative. What matters is not the absolute sign of QLIKE, but the ordering. lower is better.
We can see that the best model is not universal. Short horizons may prefer models that react quickly to shocks, while longer horizons may prefer smoother or asymmetric specifications. This is the reason that we shouldn’t pick one universal best model and use it for all the dates and we should use the best model based on recent performance for each horizon.
8) Rolling model selection
The selected forecast at date \(t\) is based only on forecast losses whose realized windows were already complete before date \(t\). We build a rolling model selector. This selector is a walk forward model selection rule.
For a forecast made on date \(s\) with horizon \(h\), the realized target is only available on
\[
s+h.
\]
So we first attach a realized availability date:
\[
a(s,h)=s+h.
\]
At current signal date \(t\), the selector is only allowed to use historical losses satisfying
\[
a(s,h)<t.
\]
For each horizon, each model has an average recent QLIKE over a rolling lookback window:
def add_realized_availability_dates(forecast_panel, market_dates): dates = pd.DatetimeIndex(pd.to_datetime(market_dates)).normalize() date_to_pos = {pd.Timestamp(d): i for i, d inenumerate(dates)} out = forecast_panel.copy() available = []for _, row in out.iterrows(): d = pd.Timestamp(row["date"]).normalize() h =int(row["horizon"]) pos = date_to_pos.get(d)if pos isNoneor pos + h >=len(dates): available.append(pd.NaT)else: available.append(pd.Timestamp(dates[pos + h]).normalize()) out["realized_available_date"] = availablereturn outdef select_model_by_rolling_past_loss(df_fc_signal, horizons=horizons, lookback=126, min_obs=40, fallback="median_ensemble", eps=eps): data = df_fc_signal.copy() data["date"] = pd.to_datetime(data["date"], errors="coerce").dt.normalize() data["horizon"] = data["horizon"].astype(int) data = add_realized_availability_dates(data, df_market["trade_date"]) data["loss_value"] = qlike_array(data["realized_var_sum"], data["forecast_var_sum"], eps=eps) rows = []for (date, horizon), current in data.groupby(["date", "horizon"], sort=True):ifint(horizon) notin horizons:continue current = current.dropna(subset=["forecast_var_sum"]).copy()if current.empty:continue past = data[ (data["horizon"] ==int(horizon))& (data["realized_available_date"].notna())& (data["realized_available_date"] < pd.Timestamp(date))& (data["loss_value"].notna())].copy()if lookback and lookback >0: keep_dates = pd.Index(sorted(past["date"].unique()))[-int(lookback):] past = past[past["date"].isin(keep_dates)] recent = ( past.groupby("model") .agg(recent_loss=("loss_value", "mean"), n_loss=("loss_value", "count")) .query("n_loss >= @min_obs") .sort_values(["recent_loss", "n_loss"], ascending=[True, False])) current_by_model = current.set_index("model", drop=False) selected_model =None selection_reason ="rolling_past_qlike" recent_loss = np.nan n_loss =0ifnot recent.empty: valid_ranked = [m for m in recent.index if m in current_by_model.index]if valid_ranked: selected_model =str(valid_ranked[0]) selected = current_by_model.loc[selected_model]ifisinstance(selected, pd.DataFrame): selected = selected.iloc[0] recent_loss =float(recent.loc[selected_model, "recent_loss"]) n_loss =int(recent.loc[selected_model, "n_loss"])if selected_model isNone: selected_model = fallback selection_reason ="fallback_median_ensemble" selected = current.iloc[0].copy() selected["forecast_var_sum"] =float(current["forecast_var_sum"].median()) selected["forecast_var_ann"] = annualization /int(horizon) *float(selected["forecast_var_sum"]) selected["forecast_vol_ann"] =float(np.sqrt(max(selected["forecast_var_ann"], 0.0))) rows.append({"date": pd.Timestamp(date),"horizon": int(horizon),"selected_model": selected_model,"selection_reason": selection_reason,"forecast_var_sum": float(selected["forecast_var_sum"]),"forecast_var_ann": float(selected["forecast_var_ann"]),"forecast_vol_ann": float(selected["forecast_vol_ann"]),"realized_var_sum": float(selected["realized_var_sum"]) if pd.notna(selected["realized_var_sum"]) else np.nan,"realized_var_ann": float(selected["realized_var_ann"]) if pd.notna(selected["realized_var_ann"]) else np.nan,"realized_vol_ann": float(selected["realized_vol_ann"]) if pd.notna(selected["realized_vol_ann"]) else np.nan,"recent_loss": recent_loss,"n_loss": n_loss, })return pd.DataFrame(rows).sort_values(["date", "horizon"]).reset_index(drop=True)selected_fc = select_model_by_rolling_past_loss(df_fc_signal, horizons=horizons, lookback=126, min_obs=40, fallback="median_ensemble")selected_model_counts = selected_fc.groupby(["horizon", "selected_model"]).size().rename("n").reset_index()display(selected_model_counts.pivot_table(index="selected_model", columns="horizon", values="n", fill_value=0))fig, ax = plt.subplots(figsize=(10.5, 4.8))plot_counts = selected_model_counts.pivot(index="selected_model", columns="horizon", values="n").fillna(0)bottom = np.zeros(len(plot_counts.columns))for model, row in plot_counts.iterrows(): vals = row.to_numpy(dtype=float) ax.bar([str(h) for h in plot_counts.columns], vals, bottom=bottom, label=model) bottom += valsax.set_title("Selected model counts by horizon")ax.set_xlabel("Forecast horizon")ax.set_ylabel("Count")ax.legend(loc="center left", bbox_to_anchor=(1.02, 0.5))plt.tight_layout()plt.show()
horizon
1
5
10
21
42
63
selected_model
egarch11_normal
76.000000
150.000000
171.000000
194.000000
203.000000
169.000000
egarch11_student
110.000000
41.000000
28.000000
1.000000
0.000000
0.000000
garch11_student
16.000000
2.000000
0.000000
0.000000
0.000000
21.000000
garch22_student
35.000000
0.000000
0.000000
0.000000
0.000000
0.000000
gjr11_normal
6.000000
77.000000
54.000000
62.000000
43.000000
110.000000
gjr11_student
212.000000
181.000000
54.000000
5.000000
0.000000
0.000000
har_rv
0.000000
0.000000
139.000000
173.000000
168.000000
94.000000
median_ensemble
41.000000
45.000000
50.000000
61.000000
82.000000
102.000000
We can see how model preference changes by horizon. EGARCH normal dominates many of the 10, 21, and 42 day selections, HAR-RV becomes important at medium and longer horizons, and the median ensemble appears frequently in periods where the rolling loss history is not yet deep enough. The one day horizon is more unstable and uses a wider mixture of models, especially GJR Student and EGARCH Student.
9) Term structure of selected forecasts and Implied volatility
Now we have two term structures that we built for SPX in 2022-2023:
the option market term structure which is computed from solving option prices for volatility for each maturity of that date
the model forecast term structure which is computed from selected realized variance forecasts across horizons for each date.
The option side is observed at maturity \(\tau\) or DTE. The model side is produced at discrete trading day horizons \(h\in\{1,5,10,21,42,63\}\). Therefore we don’t have a perfect one to one comparison, but it is still informative because both curves are annualized volatilities. If:
then implied volatility is cheap relative to the model forecast.
We show this relation in four sample dates that are kind of in different volatility states
Show code
term_iv_panel = iv_panel_all.copy()term_iv_panel = term_iv_panel[ term_iv_panel["dte_trading"].between(min(horizons), max(horizons) *365.25/252.0, inclusive="both")].copy()term_iv_panel = term_iv_panel.replace([np.inf, -np.inf], np.nan).dropna(subset=["date", "dte_trading", "atm_iv_mid"])forecast_term_dates = selected_fc.groupby("date")["horizon"].nunique()iv_term_dates = term_iv_panel.groupby("date")["expiry"].nunique()eligible_term_dates = pd.Index(sorted(set(forecast_term_dates[forecast_term_dates >=len(horizons)].index) &set(iv_term_dates[iv_term_dates >=3].index)))iflen(eligible_term_dates) ==0:print("No dates have enough IV expiries and selected forecast horizons for the term-structure comparison.")else: sample_positions = np.unique(np.linspace(0, len(eligible_term_dates) -1, min(4, len(eligible_term_dates))).round().astype(int)) term_sample_dates = [pd.Timestamp(eligible_term_dates[i]).normalize() for i in sample_positions] fig, axes = plt.subplots(2, 2, figsize=(13.5, 8.4), sharey=True) axes = axes.ravel()for ax, sample_date inzip(axes, term_sample_dates): iv_curve = term_iv_panel[term_iv_panel["date"].eq(sample_date)].sort_values("dte_trading") fc_curve = selected_fc[selected_fc["date"].eq(sample_date)].sort_values("horizon") ax.plot(iv_curve["dte_trading"], iv_curve["atm_iv_mid"], marker="o", lw=1.25, label="ATM IV by expiry") ax.plot( fc_curve["horizon"], fc_curve["forecast_vol_ann"], marker="s", lw=1.25, label="selected forecast by horizon") ax.set_title(str(sample_date.date())) ax.set_xlabel("trading-day maturity / forecast horizon") ax.set_ylabel("annualized volatility") ax.legend(loc="best")for ax in axes[len(term_sample_dates):]: ax.axis("off") plt.tight_layout() plt.show() display(selected_fc[selected_fc["date"].isin(term_sample_dates)][["date", "horizon", "selected_model", "forecast_vol_ann", "realized_vol_ann"]].head(24))
date
horizon
selected_model
forecast_vol_ann
realized_vol_ann
0
2022-01-03
1
median_ensemble
0.118018
0.009998
1
2022-01-03
5
median_ensemble
0.125970
0.142585
2
2022-01-03
10
median_ensemble
0.136869
0.162306
3
2022-01-03
21
median_ensemble
0.150552
0.185920
4
2022-01-03
42
median_ensemble
0.166578
0.206689
5
2022-01-03
63
median_ensemble
0.175724
0.212171
990
2022-09-02
1
gjr11_student
0.264238
0.065140
991
2022-09-02
5
gjr11_normal
0.239860
0.191962
992
2022-09-02
10
egarch11_normal
0.216122
0.271466
993
2022-09-02
21
egarch11_normal
0.226109
0.278505
994
2022-09-02
42
egarch11_normal
0.236715
0.268215
995
2022-09-02
63
egarch11_normal
0.237688
0.264648
1980
2023-05-04
1
gjr11_student
0.156090
0.290601
1981
2023-05-04
5
gjr11_student
0.158626
0.138261
1982
2023-05-04
10
har_rv
0.131728
0.128922
1983
2023-05-04
21
har_rv
0.143254
0.131136
1984
2023-05-04
42
har_rv
0.148594
0.117405
1985
2023-05-04
63
har_rv
0.161278
0.109553
2970
2023-12-29
1
gjr11_normal
0.103927
0.090114
2971
2023-12-29
5
gjr11_student
0.105326
0.124708
2972
2023-12-29
10
gjr11_student
0.107555
0.094939
2973
2023-12-29
21
har_rv
0.102894
0.108368
2974
2023-12-29
42
har_rv
0.111801
0.117110
2975
2023-12-29
63
har_rv
0.112708
0.110592
2022-01-03: the model forecast curve is above the ATM IV curve across much of the term structure. This means the model expected more realized movement than the option market was pricing.
2022-09-02: IV and selected forecasts are much closer, and both curves sit around a high-volatility level. This is a more balanced regime.
2023-05-04: ATM IV is generally above the model forecast curve, especially at longer maturities, which means market is pricing options based on a higher volatility relative to realized volatility forecasts.
2023-12-29: both curves are low and relatively flat, showing a calmer state.
10) Variance risk premium
The variance risk premium is the bridge between theoretical forecasting and real markets. The option market quotes implied volatility, while the econometric model forecasts realized variance. To compare them, both must be placed on the same maturity scale.
The selected forecast curve should be interpolated to the option’s trading day DTE. If the option has trading day maturity \(d\), the interpolated forecasted variance sum is
The z-score measures standardized distance from a recent mean. This is sensitive to the scale of recent VRP volatility. If VRP has been very stable, a small move can produce a large z-score. If VRP has been noisy, the same raw move may look less extreme.
The rank measures empirical extremeness. It is less sensitive to outliers and easier to interpret. A rank of 0.10 means the current VRP is in the cheapest 10% of the recent window. A rank of 0.90 means it is in the richest 10%.
Using both makes the signal more robust. The rank catches relative historical extremes, while the z-score captures standardized deviations.
Show code
def rolling_past_percentile(series, window=63, min_periods=20): values = pd.Series(series, dtype=float).to_numpy() out = np.full(len(values), np.nan, dtype=float)for i inrange(len(values)): start =max(0, i -int(window)) hist = values[start:i] hist = hist[np.isfinite(hist)]iflen(hist) >= min_periods and np.isfinite(values[i]): out[i] =float(np.mean(hist <= values[i]))return pd.Series(out, index=pd.Series(series).index)def interpolate_selected_variance_to_dte(selected_fc, iv_daily, horizons=horizons, annualization=annualization): fc = selected_fc.copy() fc["date"] = pd.to_datetime(fc["date"], errors="coerce").dt.normalize() fc["horizon"] = fc["horizon"].astype(int) iv = iv_daily.copy() iv["date"] = pd.to_datetime(iv["date"], errors="coerce").dt.normalize() rows = [] by_date = {d: g.sort_values("horizon") for d, g in fc.groupby("date")}for _, opt in iv.iterrows(): date = pd.Timestamp(opt["date"]).normalize() g = by_date.get(date) dte_trading =float(opt["dte_trading"])if g isNoneornot np.isfinite(dte_trading) or dte_trading <=0:continue g = g.dropna(subset=["forecast_var_sum"]).sort_values("horizon")if g.empty:continue x = g["horizon"].to_numpy(dtype=float) forecast_sum = g["forecast_var_sum"].to_numpy(dtype=float) realized_sum = g["realized_var_sum"].to_numpy(dtype=float) target =float(np.clip(dte_trading, np.nanmin(x), np.nanmax(x))) fc_sum_dte =float(np.interp(target, x, forecast_sum)) rv_sum_dte =float(np.interp(target, x, realized_sum)) if np.isfinite(realized_sum).any() else np.nan fc_var_ann = annualization / target * fc_sum_dte rv_var_ann = annualization / target * rv_sum_dte if np.isfinite(rv_sum_dte) else np.nan nearest = g.iloc[int(np.argmin(np.abs(x - target)))] row = opt.to_dict() row.update({"target_horizon_trading_days": target,"forecast_var_sum_dte": fc_sum_dte,"forecast_var_ann": fc_var_ann,"forecast_vol_ann": float(np.sqrt(max(fc_var_ann, 0.0))),"future_realized_var_sum": rv_sum_dte,"future_realized_var_ann": rv_var_ann,"future_realized_vol_ann": float(np.sqrt(max(rv_var_ann, 0.0))) if np.isfinite(rv_var_ann) else np.nan,"selected_model": nearest["selected_model"],"selection_reason": nearest["selection_reason"]}) rows.append(row) out = pd.DataFrame(rows).sort_values("date").reset_index(drop=True)if out.empty:return out out["atm_iv_mid"] = pd.to_numeric(out["atm_iv_mid"], errors="coerce")if out["atm_iv_mid"].dropna().median() >2.0: out["atm_iv_mid"] = out["atm_iv_mid"] /100.0 out["iv_var_ann"] = out["atm_iv_mid"] **2 out["vrp_var"] = out["iv_var_ann"] - out["forecast_var_ann"] out["vol_spread"] = out["atm_iv_mid"] - out["forecast_vol_ann"] z_window =63 min_periods =20 past_mean = out["vrp_var"].shift(1).rolling(z_window, min_periods=min_periods).mean() past_std = out["vrp_var"].shift(1).rolling(z_window, min_periods=min_periods).std(ddof=0) out["vrp_z"] = (out["vrp_var"] - past_mean) / past_std.replace(0, np.nan) out["vrp_rank"] = rolling_past_percentile(out["vrp_var"], window=z_window, min_periods=min_periods).to_numpy()return outvrp_panel = interpolate_selected_variance_to_dte(selected_fc, iv_30d_daily)required_vrp_cols = ["date", "expiry", "strike", "dte_calendar", "dte_trading", "tau", "spot", "forward","atm_iv_mid", "forecast_var_ann", "forecast_vol_ann", "future_realized_var_ann","future_realized_vol_ann", "vrp_var", "vol_spread", "vrp_z", "vrp_rank", "selected_model"]display(vrp_panel[[c for c in required_vrp_cols if c in vrp_panel.columns]].head())display(vrp_panel[["vrp_var", "vol_spread", "vrp_z", "vrp_rank"]].describe())
date
expiry
strike
dte_calendar
dte_trading
tau
spot
forward
atm_iv_mid
forecast_var_ann
forecast_vol_ann
future_realized_var_ann
future_realized_vol_ann
vrp_var
vol_spread
vrp_z
vrp_rank
selected_model
0
2022-01-03
2022-02-04
4,795.000000
32.000000
22.078029
0.087611
4,795.570000
4,792.699898
0.124809
0.023162
0.152191
0.035363
0.188049
-0.007585
-0.027383
NaN
NaN
median_ensemble
1
2022-01-04
2022-02-04
4,775.000000
31.000000
21.388090
0.084873
4,793.190000
4,791.048784
0.131639
0.021292
0.145917
0.042105
0.205194
-0.003963
-0.014278
NaN
NaN
median_ensemble
2
2022-01-05
2022-02-04
4,705.000000
30.000000
20.698152
0.082136
4,700.640000
4,699.800813
0.156209
0.035953
0.189612
0.037389
0.193364
-0.011551
-0.033403
NaN
NaN
median_ensemble
3
2022-01-06
2022-02-04
4,685.000000
29.000000
20.008214
0.079398
4,696.250000
4,695.048096
0.161918
0.030902
0.175790
0.037459
0.193542
-0.004684
-0.013871
NaN
NaN
median_ensemble
4
2022-01-07
2022-02-07
4,670.000000
31.000000
21.388090
0.084873
4,676.410000
4,674.648924
0.152411
0.027443
0.165658
0.038811
0.197005
-0.004214
-0.013248
NaN
NaN
median_ensemble
vrp_var
vol_spread
vrp_z
vrp_rank
count
496.000000
496.000000
476.000000
476.000000
mean
-0.005757
-0.012284
0.011296
0.500693
std
0.018905
0.033660
1.230328
0.310434
min
-0.164693
-0.201553
-8.985464
0.000000
25%
-0.010363
-0.028639
-0.617243
0.206349
50%
-0.003663
-0.011763
0.043781
0.507937
75%
0.002919
0.007726
0.655698
0.780303
max
0.040583
0.092577
3.523038
1.000000
We can see that the sample average VRP is negative, about -0.0057 in variance units, with a median of about -0.0036. The minimum is very negative, about -0.164, while the maximum is about 0.0406. This asymmetry means the largest dislocations in this specific sample are cheap IV states where model forecast variance exceeded implied variance by a large amount.
The rolling rank has a mean close to 0.50. This makes it useful for trading thresholds. a low rank is cheap IV, a high rank is rich IV.
VRP distribution is not symmetric. There are much larger negative spikes than positive spikes in this SPX 2 year sample. This is important because strategies that will be built with this signal may not have symmetric frequency or payoff shape. Volatility spread is easier to interpret visually but less theoretically clean. So we use VRP as the main signal.
10.2 Predicting future RV with IV and Forecasts.
We use a simple OLS model to see which of these has better prediction of future realized volatility. The regression specifications are
We use HAC standard errors because future realized variance windows overlap. If the option on date \(t\) is roughly 30 days to expiry, then the realized variance target for date \(t\) and the target for date \(t+1\) share most of the same return observations. Then ordinary errors would be too optimistic.
Using both increases \(R^2\) only slightly to about 0.520.
In the joint regression, the IV coefficient remains strongly significant, while the forecast-variance coefficient becomes statistically insignificant.
This doesn’t mean forecast is udeless. It means that, in this sample, the ATM IV level already contains most of the linear information about future realized variance. The model forecast still matters because the trading signal is relative value:
ATM IV, selected forecast volatility, and future realized volatility move together at a broad level, but they separate during stressed periods. We can also see the constructed VRP over time with negative spikes in 2022 where the model forecast was much higher than IV.
A useful way to read the signal is:
cheap IV state: \(\text{rank}_t\leq0.30 \quad\text{or}\quad z_t\leq -0.50\)
rich IV state: \(\text{rank}_t\geq0.70 \quad\text{or}\quad z_t\geq 0.50\)
11) Straddle strategy for option trading
A long straddle buys a call and a put with the same strike \(K\) and the same expiry. If the call premium is \(C\) and the put premium is \(P\), the terminal payoff is
For a call and a put with both strikes \(K=100\) and call premium \(C=6\), and put premium \(P=5\), the total premium is 11. The break even prices are therefore 89 and 111. The payoff is V shaped because the position benefits from a large move in either direction.
This payoff can be a good strategy to use VRP as a signal
When implied volatility is cheap relative to forecast volatility, the premium paid for the V-shaped payoff may be too low. That motivates a long straddle.
When implied volatility is rich relative to forecast volatility, the market may be overpaying for convexity. That motivates a short straddle.
A two-sided strategy can buy convexity in cheap-IV states and sell volatility in rich-IV states at the same time, catching all the signals made by VRP.
Show code
from quantfinlab.plotting.diagrams import plot_straddle_payofffig, ax = plot_straddle_payoff( spot=100.0, strike=100.0, call_premium=6.0, put_premium=5.0, side="long", title="Long straddle payoff: buy call + buy put at the same strike")plt.show()
The long straddle enters at the ask and exits at the bid:
Here \(u\) is the number of contracts and \(M\) is the contract multiplier. This bid/ask convention is important because it prevents the backtest from assuming that trades can be executed at mid prices. The spread cost is approximately
Warning: at least one VRP strategy has fewer than 20 trades under the default integer-contract caps.
Running documented sensitivity: cheap rank <= 0.40, rich rank >= 0.60, with larger premium/margin caps.
We first try the default thresholds and integer contract sizing. Since at least one VRP strategy produces fewer than 20 trades, we switch to looser thresholds:
AlwaysLongStraddle loses about $14,185 with a negative Sharpe around -0.40. This makes sense with the idea that systematically buying convexity is expensive when option premiums are not selectively cheap.
LongCheapConvexity earns about $6,990, but it has only 18 trades, so the result shouldn’t be trusted that much.
ShortRichVolatility earns about $33,560, has a hit rate above 73%, and reports the best Sharpe, around 1.67.
TwoSidedVRPSwitcher earns about $24,520, with a Sharpe around 0.97.
We can see that the VRP signal appears to separate states better than the unconditional long straddle benchmark. The unconditional buyer pays for convexity too often, while the VRP aware strategies trade only when the model implied relative value is worth trading on. is a controlled research overlay, not a production-ready options portfolio.
The short-rich and long convexity results should not be trusted that much because we only have a 2 year window sample and we can only say that in this sample and the regimes of this 2 years, Short volatility models would perform better. This sample may not contain the full tail risk distribution.
performance across regimes
If the VRP signal had no meaning, we might see long straddles perform randomly across cheap, neutral, and rich VRP. but the unconditional long straddle makes money in cheap VRP and loses heavily in rich VRP. This is exactly the pattern we expect and shows that a strategy without considering VRP as a signal can lose more money than a stragtegy that does.
We consider these thresholds for regime definition:
The decomposition strongly supports the interpretation of the signal:
Always Long Straddle makes money in cheap VRP states, about +$8,731, but loses heavily in rich VRP states, about -$27,110. Buying options is attractive when IV is cheap and unattractive when IV is rich.
Long Cheap Convexity earns most of its positive P&L in cheap VRP states, about +$10,060, while losing in the few neutral trades.
Short Rich Volatility earns most of its P&L in rich VRP states, about +$23,440, with additional gains in neutral states.
Two Sided VRP Switcher is positive across all the regimes, but its largest contribution comes from rich VRP states, which means the short volatility side is doing much of the work.
Show code
ifnot all_equity.empty: fig, axes = plt.subplots(3, 2, figsize=(14.5, 12.0)) axes = axes.ravel()for name, g in all_equity.groupby("strategy"): axes[0].plot(g["date"], g["nav"], lw=1.1, label=name) axes[1].plot(g["date"], g["drawdown"], lw=1.0, label=name) axes[0].set_title("Strategy NAV comparison") axes[1].set_title("Strategy drawdown comparison") axes[0].legend(loc="best") axes[1].legend(loc="best") axes[2].plot(vrp_panel["date"], vrp_panel["vrp_rank"], lw=0.9, label="VRP rank") axes[2].axhline(0.30, ls="--", lw=0.8, c="tab:green") axes[2].axhline(0.70, ls="--", lw=0.8, c="tab:red")ifnot all_trades.empty: marker_dates = all_trades["entry_date"].dropna().unique() marker_panel = vrp_panel[vrp_panel["date"].isin(marker_dates)] axes[2].scatter(marker_panel["date"], marker_panel["vrp_rank"], s=18, c="black", alpha=0.55, label="entries") axes[2].set_title("VRP rank with trade-entry markers") axes[2].legend(loc="best")ifnot all_trades.empty:for name, g in all_trades.groupby("strategy"): axes[3].hist(g["net_pnl"], bins=30, alpha=0.45, label=name) axes[3].set_title("Trade P&L distribution by strategy") axes[3].legend(loc="best") decile_plot = all_trades.dropna(subset=["vrp_decile"]).groupby(["vrp_decile"], observed=True)["net_pnl"].sum() axes[4].bar([str(x) for x in decile_plot.index], decile_plot.to_numpy(dtype=float)) axes[4].tick_params(axis="x", rotation=45) axes[4].set_title("Total P&L by VRP rank decile") switcher = all_trades[all_trades["strategy"] =="TwoSidedVRPSwitcher"]ifnot switcher.empty: side_pnl = switcher.groupby("side")["net_pnl"].sum() axes[5].bar(side_pnl.index.astype(str), side_pnl.to_numpy(dtype=float)) axes[5].set_title("TwoSidedVRPSwitcher P&L by side")else: axes[5].axis("off")else:for ax in axes[3:]: ax.axis("off") plt.tight_layout() plt.show()
Implementation of the project on BTC-Deribit options chain with QuantFinLab
We use an option chain of Bitcoin for the first 9 months of 2023 to show the implementation with our library. We picked it because it’s european and doesn’t include dividends and can be closer to the implementation on SPX. But BTC requires several adjustments.
First, crypto trades every calendar day, so the annualization changes from 252 to 365
Second, BTC option quotes may be represented in different units, so we detect and converts option prices into USD equivalent before building IV.
Third, for this case we use fractional units and a contract multiplier of 1 for the backtest, because the BTC option scale is different from SPX index options.
The BTC ATM-IV cleaning shows we started with about 258,234 raw rows and ended up with 1,228 final ATM rows. The largest filters are the DTE and moneyness filters, which is expected for an option chain dataset with many strikes and expiries.
The BTC forecast ranking is different from SPX. based on QLIKE the overal first is GJR normal first, followed by EGARCH normal and GARCH normal. At the 21-day horizon, GJR normal and EGARCH normal are again the best by QLIKE. The Student-\(t\) models perform worse in this BTC experiment. Parameter instability, small sample length, crypto regime shifts, and noisy option implied inputs can all change the ranking.
Based on The selected model count table, EGARCH normal is heavily used at one-day horizons, while HAR-RV dominates the 63-day horizon. Short horizon forecasts often need fast reaction to recent shocks, while longer horizons benefit from smoother realized volatility memory.
The BTC VRP and strategy tables are less trustable than the SPX results because the sample is shorter and the option microstructure is different.
Always Long Straddle loses across the scenarios, with about -0.26% return on initial NAV in the displayed loose/base/strict rows.
Long Cheap Convexity is also negative in the scenarios.
Short Rich Volatility is slightly positive in loose/base scenarios but often has very few trades.
Two Sided VRP Switcher is negative in the scenario set, although the side breakdown shows that short trades can still be profitable.
