HAR-Volatilitätsprognose und Rauigkeitsmessung

# HAR and Rough Volatility


# HAR and Rough Volatility

**Chapter 9 | Section 9.3**

**Docker image**: `ml4t`

The GARCH model in the previous notebook has one memory: a single decay rate that
governs how long a shock lasts. This notebook covers two things that memory cannot
express. **HAR** replaces the single decay with three horizons at once, on the argument
that the people trading at a daily, weekly and monthly rhythm are different people.
**Roughness** asks a different question entirely: not how long volatility remembers, but
how jagged its path is, and the answer for real markets is far jaggeder than the models
that were standard for thirty years assume.

Both need a measurement of volatility to work on, and getting that measurement right is
the first third of the notebook.

**Learning objectives**

- Compute volatility from within-session data, and use it to check what the estimators
  available from daily bars alone are actually measuring.
- Fit a model of volatility across three horizons, and read its coefficients as a
  statement about which horizon is driving the current level.
- Correct the standard errors of that fit for the overlap its own regressors contain.
- Estimate how jagged a series is over time, apply it to returns and to volatility, and
  read the two answers as the different things they are.

**Book reference**

Chapter 9, Section 9.3 (Volatility Features).

**Prerequisites**

`08_garch_volatility` for conditional volatility and for what persistence means.
Chapter 8 introduced the range-based estimators this notebook tests.

## Setup

```python
"""HAR and rough volatility - multi-horizon volatility and the Hurst exponent."""

import warnings

import exchange_calendars as xcals
import matplotlib.dates as mdates
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import statsmodels.api as sm
from arch import arch_model
from IPython.display import display
from ml4t.engineer.features.regime import hurst_exponent
from ml4t.engineer.features.volatility import (
    garman_klass_volatility,
    parkinson_volatility,
    realized_volatility,
    rogers_satchell_volatility,
    yang_zhang_volatility,
)
from pandas.errors import PerformanceWarning
from scipy import stats

from data import load_etfs, load_nasdaq100_bars
from utils.style import COLORS, FIGSIZE, show_with_alt

# Building the exchange session index applies a date offset element by element, which pandas
# reports as slow twenty times over. It is a speed notice about one calendar build.
warnings.filterwarnings(
    "ignore", category=PerformanceWarning, module="exchange_calendars.exchange_calendar"
)
```

```python
START_DATE = "2006-01-01"
END_DATE = "2024-12-31"
INTRADAY_SYMBOL = "AAPL"
DAILY_SYMBOL = "SPY"
MINIMUM_BARS_PER_SESSION = 360  # a full session is about 390 one-minute bars
EXCHANGE = "XNYS"  # the calendar that says which dates are sessions at all
SESSIONS_PER_YEAR = 252
ROLLING_WINDOW = 20  # sessions in every rolling average below: about one month
```

# Part 1: what "volatility" is being measured

Everything in this notebook models a series called realized volatility, and there are
several different things that name refers to. The definition with the strongest claim to
it is computed from within-session data: sum the squared returns of every minute in a
session, and the total estimates that session's variance with an error that falls as the
bars get finer.

$$RV_t = \sum_{i=1}^{M} r_{t,i}^2$$

It measures the **trading session** and nothing else, and that is the first thing to be
clear about, because a daily estimator's target is the whole day. This section
establishes the relation between the two on a symbol where both are available, and turns
up two data traps on the way.

```python
minute_bars = (
    load_nasdaq100_bars(symbols=[INTRADAY_SYMBOL], include_microstructure=True, lazy=True)
    .filter((pl.col("time") >= "09:30") & (pl.col("time") < "16:00"))
    .select(["date", "time", "first_trade_price", "last_trade_price"])
    .collect()
    .sort(["date", "time"])
)

print(
    f"{INTRADAY_SYMBOL} minute bars: {minute_bars.height:,} rows over "
    f"{minute_bars['date'].n_unique()} sessions, "
    f"{minute_bars['date'].min()} to {minute_bars['date'].max()}"
)
```

Two choices in the cell above are worth stating rather than leaving in the code. The
returns are taken from the **last trade price** of each minute rather than from the
minute's volume-weighted average. A volume-weighted average is an average over the
minute, so a series of them is smoother than the price actually was, and squared returns
built from it understate the variance. On this sample the difference is not small.

And the bars cover 09:30 to 16:00 only, so the sum below is the variance of the trading
session. What happens between one session's close and the next one's open is a separate
quantity, measured separately here.

Three details in the cell below are load-bearing. The session's open and close come from
the bars as loaded, not from a frame reduced to the minutes that have a previous close;
an open read from such a frame would be the second minute's price, and the overnight
return would then swallow the opening minute. The opening minute keeps its own
open-to-close return in the variance sum for the same reason: without it the sum covers
one minute less than the session return it is compared against.

And an overnight return is only an overnight return if the two rows are consecutive
sessions of the exchange. Bar counts cannot establish that, and neither can the number of
calendar days between the dates: a Monday-to-Wednesday pair is three days apart and still
contains all of Tuesday's trading if Tuesday is missing from the data. What settles it is
the exchange calendar, which says which dates are sessions independently of what the data
holds. Each pair is required to be one session apart on that calendar, which keeps
holiday weekends and drops a missing weekday. On a reduced test fixture, which keeps
separated blocks of sessions, an unchecked gap can be weeks long.

```python
boundaries = (
    minute_bars.group_by("date")
    .agg(
        session_open=pl.col("first_trade_price").first(),
        session_close=pl.col("last_trade_price").last(),
        bars=pl.len(),
    )
    .sort("date")
)

variances = (
    minute_bars.with_columns(
        # The first bar of a session has no previous close, and dropping it leaves the
        # opening minute out of the sum while the aggregate open-to-close return contains
        # it. Its own open-to-close return is what belongs there.
        minute_return=pl.when(pl.col("last_trade_price").log().diff().over("date").is_null())
        .then((pl.col("last_trade_price") / pl.col("first_trade_price")).log())
        .otherwise(pl.col("last_trade_price").log().diff().over("date"))
    )
    .group_by("date")
    .agg(minute_variance=pl.col("minute_return").pow(2).sum())
)

exchange_sessions = xcals.get_calendar(EXCHANGE).sessions_in_range(
    minute_bars["date"].min(), minute_bars["date"].max()
)
session_numbers = pl.DataFrame(
    {
        "date": pl.Series("date", exchange_sessions).cast(pl.Date),
        "session_number": np.arange(len(exchange_sessions), dtype=np.int64),
    }
)

intraday = (
    boundaries.join(variances, on="date", how="left")
    .join(session_numbers, on="date", how="inner")
    .sort("date")
    .with_columns(
        previous_close=pl.col("session_close").shift(1),
        previous_bars=pl.col("bars").shift(1),
        sessions_apart=pl.col("session_number") - pl.col("session_number").shift(1),
        days_apart=(pl.col("date") - pl.col("date").shift(1)).dt.total_days(),
    )
    .with_columns(
        session_return=(pl.col("session_close") / pl.col("session_open")).log(),
        overnight_return=(pl.col("session_open") / pl.col("previous_close")).log(),
        close_to_close_return=(pl.col("session_close") / pl.col("previous_close")).log(),
    )
    .filter(
        (pl.col("bars") >= MINIMUM_BARS_PER_SESSION)
        & (pl.col("previous_bars") >= MINIMUM_BARS_PER_SESSION)
        & (pl.col("sessions_apart") == 1)
    )
    .drop_nulls()
)

print(f"Sessions kept: {intraday.height:,} of {boundaries.height}")
print(f"Exchange sessions in the sample period: {len(exchange_sessions):,}")
print(f"Longest calendar gap among the kept pairs: {intraday['days_apart'].max()} days")
```

## The first trap: these prices are not adjusted

A close-to-close return computed from raw trade prices crosses every corporate action as
though it were a price move. The check below finds any session whose return is too large
to be one, and on this symbol and period it finds exactly one.

```python
IMPLAUSIBLE_RETURN = 0.25  # a one-session log return larger than this is a corporate action

implausible = intraday.filter(pl.col("close_to_close_return").abs() > IMPLAUSIBLE_RETURN)
print(f"Sessions with a return beyond {IMPLAUSIBLE_RETURN:.0%}: {implausible.height}")
for row in implausible.iter_rows(named=True):
    print(f"  {row['date']}: {row['close_to_close_return']:+.3f} in logs")

clean = intraday.filter(pl.col("close_to_close_return").abs() <= IMPLAUSIBLE_RETURN)
```

One session, and it is the four-for-one share split. The intraday variance of that
session is unaffected, because a split happens between sessions and every minute return
inside the session is a real price move; what is corrupted is anything computed across
the boundary. That asymmetry is worth carrying: the same raw data can be sound for one
statistic and unusable for another, and which is which depends on whether the statistic
crosses a session boundary.

## The second trap: a session is not a day

A close-to-close return is the overnight move plus the session that follows it, exactly,
because both are logarithms. So the identity below holds term by term with no
approximation:

$$\mathbb{E}[r_{cc}^2] = \mathbb{E}[r_{on}^2] + \mathbb{E}[r_{oc}^2]
  + 2\,\mathbb{E}[r_{on} r_{oc}]$$

where $r_{oc}$ is the **aggregate** open-to-close return of the session. That is not the
same quantity as the realized variance measured from minute bars: the realized variance
is the sum of squared minute returns and the aggregate is the square of their sum, so
they differ by every cross-product between minutes. Both are reported below, separately,
because they answer different questions and conflating them is how a decomposition comes
out looking exact when it is not.

```python
annualize = np.sqrt(SESSIONS_PER_YEAR)

overnight_variance = float((clean["overnight_return"] ** 2).mean())
session_aggregate_variance = float((clean["session_return"] ** 2).mean())
cross_product = float((clean["overnight_return"] * clean["session_return"]).mean())
close_to_close_variance = float((clean["close_to_close_return"] ** 2).mean())
realized_variance = float(clean["minute_variance"].mean())

display(
    pd.DataFrame(
        [
            {"term": "overnight, squared", "value": overnight_variance},
            {"term": "session open to close, squared", "value": session_aggregate_variance},
            {"term": "twice their cross-product", "value": 2 * cross_product},
            {
                "term": "the three added",
                "value": overnight_variance + session_aggregate_variance + 2 * cross_product,
            },
            {"term": "close to close, squared", "value": close_to_close_variance},
        ]
    )
)
print(
    f"Share of the daily square that happens overnight: {overnight_variance / close_to_close_variance:.1%}"
)
print(f"Session variance minute by minute: {np.sqrt(realized_variance) * annualize:.4f} annualized")
print(
    "Session variance from the aggregate open-to-close return: "
    f"{np.sqrt(session_aggregate_variance) * annualize:.4f} annualized"
)
print(f"Ratio of the two: {realized_variance / session_aggregate_variance:.3f}")
```

The identity closes to the last digit, which it must; what it buys is the size of each
term. The overnight square is a large fraction of the daily square for a symbol whose
market is open six and a half hours out of twenty-four, and the cross-product is small
beside the other two, so the two components are close to uncorrelated.

The last ratio is the separate question. Summing squared minute returns and squaring the
aggregate open-to-close return are two estimates of the same session's variance and they
do not agree: the minute-by-minute sum picks up movement that reverses before the close
and the aggregate does not. Neither is wrong, and which one a model wants depends on
whether reversals inside a session are risk it is exposed to. A position held through the
session sees the aggregate; one traded inside it sees the sum.

Either way, a comparison between an intraday measurement and a daily estimator has to
account for the overnight term, and a strategy holding overnight is exposed to a risk no
intraday measurement sees.

```python
fig, axes = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)

sessions = clean["date"].to_list()
ax = axes[0]
ax.plot(
    sessions,
    clean["minute_variance"].sqrt() * annualize,
    linewidth=0.7,
    color=COLORS["blue"],
    label="Trading session",
)
ax.plot(
    sessions,
    clean["overnight_return"].abs() * annualize,
    linewidth=0.5,
    alpha=0.7,
    color=COLORS["amber"],
    label="Overnight",
)
ax.set_ylabel("Annualized volatility")
ax.set_title("Two components of the same day")
ax.legend(fontsize=7)

ax = axes[1]
ax.plot(
    sessions,
    (clean["minute_variance"].sqrt() * annualize).rolling_mean(ROLLING_WINDOW),
    linewidth=1,
    color=COLORS["blue"],
    label="Trading session",
)
ax.plot(
    sessions,
    (clean["overnight_return"].abs() * annualize).rolling_mean(ROLLING_WINDOW),
    linewidth=1,
    color=COLORS["amber"],
    label="Overnight",
)
ax.set_ylabel("Annualized volatility")
ax.set_xlabel("Session")
ax.set_title(f"The same two, averaged over {ROLLING_WINDOW} sessions")
ax.legend(fontsize=7)

fig.suptitle("The market is shut for most of the day, and it moves then too")
show_with_alt(
    fig,
    "Two stacked panels sharing a time axis over two years. The top draws the session "
    "volatility measured from minute bars against the absolute overnight move, both noisy "
    "day to day, with the session series generally the larger and both spiking in March "
    "2020. The bottom draws the same two averaged over twenty sessions, where the session "
    "series stays above the overnight one throughout and the two move together.",
)
```

## Which daily estimator to use

Ranking estimators by accuracy needs a target measured the same way, and the section
above is the reason there is not one here: the intraday measurement covers the session
and every daily estimator covers something else.

What is measured below instead is each estimator's **dispersion around its own trailing
average**, session by session on the ETF panel's adjusted daily bars. Be clear about what
that is and is not. It is not the estimator's sampling error: the underlying volatility
genuinely moves within twenty sessions, so real movement is in the number, and the
trailing average in the denominator carries the estimator's own noise as well. Nor are
the four measuring the same interval, since the close-to-close estimator spans a session
boundary and the other three do not.

What it does show is how jagged each column would be if a model read it, which is a
property of the column rather than of the estimator, and is the thing that decides
whether a feature needs smoothing before use.

```python
def range_volatility(frame: pl.DataFrame) -> pl.DataFrame:
    """Per-session Garman-Klass volatility, annualized, plus the close-to-close return."""
    return frame.with_columns(
        log_return=pl.col("close").log().diff(),
        volatility=(
            (
                0.5 * (pl.col("high") / pl.col("low")).log().pow(2)
                - (2 * np.log(2) - 1) * (pl.col("close") / pl.col("open")).log().pow(2)
            ).clip(lower_bound=0)
            * SESSIONS_PER_YEAR
        ).sqrt(),
    ).drop_nulls()


spy = (
    load_etfs(symbols=[DAILY_SYMBOL])
    .filter(pl.col("timestamp") >= pl.lit(START_DATE).str.to_date())
    .filter(pl.col("timestamp") <= pl.lit(END_DATE).str.to_date())
    .sort("timestamp")
    .select(["timestamp", "open", "high", "low", "close"])
)

per_session = spy.with_columns(
    close_to_close=(pl.col("close").log().diff().pow(2) * SESSIONS_PER_YEAR).sqrt(),
    parkinson=(
        (pl.col("high") / pl.col("low")).log().pow(2) / (4 * np.log(2)) * SESSIONS_PER_YEAR
    ).sqrt(),
    garman_klass=(
        (
            0.5 * (pl.col("high") / pl.col("low")).log().pow(2)
            - (2 * np.log(2) - 1) * (pl.col("close") / pl.col("open")).log().pow(2)
        ).clip(lower_bound=0)
        * SESSIONS_PER_YEAR
    ).sqrt(),
    rogers_satchell=(
        (
            (pl.col("high") / pl.col("close")).log() * (pl.col("high") / pl.col("open")).log()
            + (pl.col("low") / pl.col("close")).log() * (pl.col("low") / pl.col("open")).log()
        ).clip(lower_bound=0)
        * SESSIONS_PER_YEAR
    ).sqrt(),
).drop_nulls()

ESTIMATORS = ["close_to_close", "parkinson", "garman_klass", "rogers_satchell"]

efficiency = pd.DataFrame(
    [
        {
            "estimator": name,
            "mean, annualized": per_session[name].mean(),
            "dispersion around its trailing average": float(
                (per_session[name] / per_session[name].rolling_mean(ROLLING_WINDOW))
                .drop_nulls()
                .std()
            ),
            "sessions estimated at zero": int((per_session[name] < 1e-6).sum()),
        }
        for name in ESTIMATORS
    ]
).sort_values("dispersion around its trailing average")
display(efficiency)
```

```python
fig, ax = plt.subplots(figsize=FIGSIZE["single"])
smooth = per_session.select(
    [pl.col("timestamp")]
    + [pl.col(name).rolling_mean(ROLLING_WINDOW).alias(name) for name in ESTIMATORS]
).drop_nulls()
for name, color in zip(
    ESTIMATORS, [COLORS["recede"], COLORS["amber"], COLORS["blue"], COLORS["copper"]]
):
    ax.plot(smooth["timestamp"], smooth[name], linewidth=0.8, alpha=0.9, color=color, label=name)
ax.set_ylabel("Annualized volatility")
ax.set_xlabel("Session")
ax.set_title(f"Four estimators of the same thing, averaged over {ROLLING_WINDOW} sessions")
ax.legend(fontsize=7)
show_with_alt(
    fig,
    "Four volatility estimators for the same symbol, each averaged over twenty sessions and "
    "drawn over the whole sample. All four rise and fall together at the same dates and are "
    "sit almost on top of one another, with the close-to-close series the noisiest and "
    "reaching the highest peaks and the three range-based ones tighter and a little "
    "below it there.",
)
```

The close-to-close column is the one that separates from the rest, on both of the
quantities reported. It is the most dispersed around its own average, and it is the only
one that puts a non-trivial number of sessions at zero: a session that closed where it
opened has no close-to-close movement to read, whatever happened in between. Reading the
high and the low is what removes that failure, and Chapter 8 is where the efficiency of
each formula is derived rather than measured.

The mean column ranks them on level, and there the ordering says nothing about which is
right. Each formula assumes something about the process, they cover different parts of
the day, and this notebook has no target to check any of them against. Read the level as
a scale each estimator carries, and never mix two of them in one feature without
rescaling.

## The library's versions, and the one that sees the gap

`ml4t.engineer.features.volatility` supplies all five as Polars expressions. They take a
window rather than working session by session, and they average variances before taking
the square root where the cells above take the square root first, so the two do not agree
exactly even where they implement the same formula. The difference is the usual one
between a mean of roots and a root of means, and it is small next to the difference
between estimators.

The fifth estimator is the reason to run this: **Yang-Zhang** adds an overnight component
to an intraday one, so it is the only one of the five whose target is the whole day that
Part 1 decomposed. Its level should sit above the others for exactly that reason.

```python
library = (
    per_session.with_columns(log_return=pl.col("close").log().diff())
    .with_columns(
        library_close=realized_volatility("log_return", period=ROLLING_WINDOW),
        library_parkinson=parkinson_volatility("high", "low", period=ROLLING_WINDOW),
        library_garman_klass=garman_klass_volatility(
            "open", "high", "low", "close", period=ROLLING_WINDOW
        ),
        library_rogers_satchell=rogers_satchell_volatility(
            "open", "high", "low", "close", period=ROLLING_WINDOW
        ),
        library_yang_zhang=yang_zhang_volatility(
            "open", "high", "low", "close", period=ROLLING_WINDOW
        ),
    )
    .drop_nulls()
)

display(
    pd.DataFrame(
        [
            {
                "expression": name.removeprefix("library_"),
                "mean, annualized": library[name].mean(),
                "sees the overnight gap": name.endswith(("close", "yang_zhang")),
            }
            for name in library.columns
            if name.startswith("library_")
        ]
    )
)
```

Yang-Zhang and the close-to-close estimator are the two that cross a session boundary,
and they are the two whose target is the day rather than the session. The three purely
intraday estimators are measuring something smaller and should not be compared against
them on level, which is the same point Part 1 made with the decomposition and is the
reason a feature set should not mix them.

# Part 2: three horizons instead of one

The **heterogeneous autoregressive** model (Corsi, 2009) predicts tomorrow's volatility
from three averages of today's: yesterday's, the last week's, and the last month's.

$$RV_{t+1} = c + \beta_d\,RV^{(d)}_t + \beta_w\,RV^{(w)}_t + \beta_m\,RV^{(m)}_t + \varepsilon_{t+1}$$

The argument for that shape is economic rather than statistical. Traders operating at
different frequencies watch volatility at different horizons and act on what they watch,
so the volatility a market produces is a superposition of their responses rather than a
single decay. The model is a linear regression, which makes it far easier to fit than
GARCH and far easier to read: each coefficient is the weight one horizon carries.

The regressors are built from a range-based estimator for the reason Part 1 established,
and each is lagged so that nothing on the right-hand side of the regression is known
later than the session it predicts from.

```python
HORIZONS = {"daily": 1, "weekly": 5, "monthly": 22}

spy = (
    load_etfs(symbols=[DAILY_SYMBOL])
    .filter(pl.col("timestamp") >= pl.lit(START_DATE).str.to_date())
    .filter(pl.col("timestamp") <= pl.lit(END_DATE).str.to_date())
    .sort("timestamp")
    .select(["timestamp", "open", "high", "low", "close"])
)


def range_volatility(frame: pl.DataFrame) -> pl.DataFrame:
    """Per-session Garman-Klass volatility, annualized, plus the close-to-close return."""
    return frame.with_columns(
        log_return=pl.col("close").log().diff(),
        volatility=(
            (
                0.5 * (pl.col("high") / pl.col("low")).log().pow(2)
                - (2 * np.log(2) - 1) * (pl.col("close") / pl.col("open")).log().pow(2)
            ).clip(lower_bound=0)
            * SESSIONS_PER_YEAR
        ).sqrt(),
    ).drop_nulls()


def har_frame(frame: pl.DataFrame) -> pl.DataFrame:
    """The three lagged averages and the next session's volatility, as one table."""
    return (
        frame.with_columns(
            [
                pl.col("volatility").shift(1).rolling_mean(length).alias(name)
                for name, length in HORIZONS.items()
            ]
        )
        .rename({"volatility": "target"})
        .drop_nulls()
    )


spy_volatility = range_volatility(spy)
har = har_frame(spy_volatility)
print(f"{DAILY_SYMBOL}: {spy_volatility.height:,} sessions, {har.height:,} usable for the model")
display(har.select(["timestamp", *HORIZONS, "target"]).tail(3))
```

## Fitting it, and correcting the standard errors

The regression is ordinary least squares. Its textbook standard errors assume the
residuals are homoscedastic and serially uncorrelated, and neither assumption is safe on
a volatility series: volatility is heteroscedastic almost by definition, and a
misspecified dynamic leaves autocorrelation behind in the residuals.

Note what the overlapping regressors do and do not imply. The weekly and monthly
regressors are rolling averages, so consecutive rows share most of their inputs; that
makes the *regressors* correlated across rows, which ordinary least squares handles.
It does not by itself make the *residuals* correlated: a correctly specified model of
this shape can have independent innovations. So overlap is a reason to check rather than
a proof of a problem.

Four standard errors for the same coefficients are computed below, and the reason there
are four rather than two is that isolating one relaxation requires holding the other
estimator choices fixed. The textbook ones assume homoscedastic, serially uncorrelated
residuals. `HC3` relaxes the first assumption with a leverage adjustment. The last two
are the same **Newey-West** estimator run twice, once with a zero-lag window and once
over a window set to the longest horizon in the model, which is where any dependence the
regressors' overlap could induce would reach. Those two differ in the lag window and in
nothing else, so their ratio is what the lagged covariance terms contribute; `HC3`
against the textbook column is a separate reading of heteroscedasticity alone.

```python
NEWEY_WEST_LAGS = max(HORIZONS.values())


def fit_har(frame: pl.DataFrame):
    """OLS with Newey-West standard errors over the longest horizon's overlap."""
    design = sm.add_constant(frame.select(list(HORIZONS)).to_pandas())
    return sm.OLS(frame["target"].to_numpy(), design).fit(
        cov_type="HAC", cov_kwds={"maxlags": NEWEY_WEST_LAGS}
    )


design = sm.add_constant(har.select(list(HORIZONS)).to_pandas())
target = har["target"].to_numpy()

har_fit = fit_har(har)
plain_fit = sm.OLS(target, design).fit()
leverage_fit = sm.OLS(target, design).fit(cov_type="HC3")
zero_lag_fit = sm.OLS(target, design).fit(cov_type="HAC", cov_kwds={"maxlags": 0})

display(
    pd.DataFrame(
        {
            "coefficient": har_fit.params,
            "standard error, textbook": plain_fit.bse,
            "standard error, HC3": leverage_fit.bse,
            "standard error, Newey-West at zero lags": zero_lag_fit.bse,
            f"standard error, Newey-West at {NEWEY_WEST_LAGS} lags": har_fit.bse,
            "ratio, HC3 to textbook": leverage_fit.bse / plain_fit.bse,
            "ratio, the two Newey-West columns": har_fit.bse / zero_lag_fit.bse,
        }
    )
)
print(f"R-squared: {har_fit.rsquared:.4f}")
```

Nothing about the coefficients changed. What changed is how much confidence the fit
reports in them, and the two ratio columns say where that change comes from.

The first ratio is heteroscedasticity, read on its own. The second is the lag window, read
with everything else held fixed. A second ratio near one says the lagged covariance terms
add little to the zero-lag version; it is a statement about the estimate, not about the
residuals, because what Newey-West sums over the window is the covariance of the
regressor-residual products rather than of the residuals themselves. Serially correlated
residuals whose products with these regressors happen to be close to uncorrelated would
produce the same ratio.

The coefficients themselves are the reason to fit HAR rather than GARCH. Each one is the
weight the model puts on one horizon, so their relative sizes say which horizon is
carrying the current level, and that is a statement a person can act on. GARCH answers
the same forecasting question with one decay rate and no such decomposition.

Two things to watch when reading them. They are not constrained to be positive, and a
negative coefficient on one horizon usually means it is collinear with another rather
than that volatility at that horizon predicts less volatility next session. And they sum
to close to one on a persistent series, so a rise in one is generally a fall in another;
read the three together rather than one at a time.

## Against GARCH, out of sample

Both models are fitted on the same training block and asked for one-step-ahead forecasts
over the same test block. They are not forecasting exactly the same quantity, which is
the first thing the comparison has to state: HAR targets the range-based volatility it
was built from, and GARCH targets the variance of close-to-close returns. The levels
therefore differ for a reason that has nothing to do with forecast quality, and the error
measures below carry that difference.

```python
TRAIN_FRACTION = 0.7

split = int(har.height * TRAIN_FRACTION)
train, test = har.head(split), har.tail(har.height - split)

har_out_of_sample = fit_har(train).predict(sm.add_constant(test.select(list(HORIZONS)).to_pandas()))

returns_percent = pd.Series(
    spy_volatility["log_return"].to_numpy() * 100,
    index=pd.DatetimeIndex(spy_volatility["timestamp"].to_list()),
)
split_timestamp = pd.Timestamp(test["timestamp"][0])
garch_fit = arch_model(returns_percent, mean="Constant", vol="GARCH", p=1, q=1).fit(
    disp="off", last_obs=split_timestamp
)
origin = returns_percent.index[returns_percent.index.get_loc(split_timestamp) - 1]
garch_out_of_sample = (
    np.sqrt(garch_fit.forecast(horizon=1, start=origin, reindex=False).variance["h.1"].to_numpy())
    * np.sqrt(SESSIONS_PER_YEAR)
    / 100
)[: test.height]

actual = test["target"].to_numpy()
display(
    pd.DataFrame(
        [
            {
                "model": name,
                "root mean squared error": float(np.sqrt(np.mean((forecast - actual) ** 2))),
                "mean absolute error": float(np.mean(np.abs(forecast - actual))),
                "mean forecast": float(np.mean(forecast)),
                "correlation with the target": float(np.corrcoef(forecast, actual)[0, 1]),
            }
            for name, forecast in [
                ("HAR", har_out_of_sample.to_numpy()),
                ("GARCH", garch_out_of_sample),
            ]
        ]
        + [{"model": "the target itself", "mean forecast": float(np.mean(actual))}]
    )
)
```

```python
fig, axes = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)

test_sessions = test["timestamp"].to_list()
ax = axes[0]
ax.plot(test_sessions, actual, linewidth=0.5, alpha=0.6, color=COLORS["neutral"], label="Target")
ax.plot(test_sessions, har_out_of_sample, linewidth=0.9, color=COLORS["blue"], label="HAR")
ax.plot(test_sessions, garch_out_of_sample, linewidth=0.9, color=COLORS["amber"], label="GARCH")
ax.set_ylabel("Annualized volatility")
ax.set_title("One-step-ahead forecasts, parameters fitted before the block")
ax.legend(fontsize=7)

ax = axes[1]
ax.plot(test_sessions, har_out_of_sample - actual, linewidth=0.5, color=COLORS["blue"], label="HAR")
ax.plot(
    test_sessions, garch_out_of_sample - actual, linewidth=0.5, color=COLORS["amber"], label="GARCH"
)
ax.axhline(0, color=COLORS["neutral"], linestyle="--", linewidth=0.5)
ax.set_ylabel("Forecast minus target")
ax.set_xlabel("Session")
ax.set_title("Errors, on the same scale")
ax.legend(fontsize=7)

fig.suptitle("Two models of the same thing, forecasting two slightly different things")
show_with_alt(
    fig,
    "Two stacked panels over the test block. The top draws the target volatility against "
    "the HAR and GARCH forecasts: both track its large movements, the HAR line sitting "
    "closer to the target and the GARCH line at a visibly higher level throughout. The "
    "bottom draws each forecast minus the target, with the GARCH errors offset above zero "
    "and both widening at the same dates.",
)
```

Read the mean forecast column before the error columns. The two models sit at different
levels because they are estimating different quantities, and a squared error between a
forecast and a target measured a different way charges the model for that difference as
though it were inaccuracy. The correlation column is the part of the comparison that
the level difference does not touch, because a correlation is invariant to it.

So the comparison establishes less than it looks like it does. What it does establish is
that a linear regression on three lagged averages tracks realized volatility about as
well as a fitted GARCH does, at a fraction of the machinery, which is the reason HAR is
the standard baseline in the realized-volatility literature.

## The coefficients as features

Refitting on a moving window turns the three weights into three slowly varying columns,
and their movement is a statement about which horizon is driving volatility now. The
window is long because the coefficients are estimated rather than observed, and the step
is a month because refitting daily would produce three columns whose variation is mostly
estimation noise.

```python
REFIT_WINDOW = 504  # sessions each fit reads: about two years
REFIT_STEP = 22  # sessions between fits: about one month

rolling_rows = []
for start in range(0, har.height - REFIT_WINDOW, REFIT_STEP):
    window = har.slice(start, REFIT_WINDOW)
    fit = fit_har(window)
    rolling_rows.append(
        {
            "timestamp": window["timestamp"][-1],
            **{name: fit.params[name] for name in HORIZONS},
            "r_squared": fit.rsquared,
        }
    )

rolling_har = pd.DataFrame(rolling_rows).set_index("timestamp")
print(f"Refits: {len(rolling_har)}, each on {REFIT_WINDOW} sessions")
```

```python
fig, axes = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)

ax = axes[0]
for name, color in zip(HORIZONS, [COLORS["blue"], COLORS["amber"], COLORS["copper"]]):
    ax.plot(rolling_har.index, rolling_har[name], linewidth=1, color=color, label=name)
ax.axhline(0, color=COLORS["neutral"], linestyle="--", linewidth=0.5)
ax.set_ylabel("Coefficient")
ax.set_title("Which horizon carries the weight, and when")
ax.legend(fontsize=7)

ax = axes[1]
ax.fill_between(rolling_har.index, 0, rolling_har["r_squared"], alpha=0.3, color=COLORS["blue"])
ax.plot(rolling_har.index, rolling_har["r_squared"], linewidth=1, color=COLORS["blue"])
ax.set_ylabel("R-squared")
ax.set_xlabel("Session the window ends")
ax.set_title("How much of the variation the fit explains")

fig.suptitle(f"Refitted every {REFIT_STEP} sessions, each fit reading {REFIT_WINDOW}")
show_with_alt(
    fig,
    "Two stacked panels over the refit dates. The top draws the three HAR coefficients, "
    "which move substantially and cross each other repeatedly; the weekly one is the "
    "largest over much of the sample, the daily one overtakes it in places, and the monthly "
    "one dips below zero in a few windows. The bottom draws the fit's R-squared, which "
    "rises and falls between about a sixth and three quarters across the sample.",
)
```

## The shape of the term structure

The ratio of the short horizon to the long one is a single number saying whether
volatility right now is above or below where it has been. Above one, the recent past has
been more volatile than the month; below one, less. It is bounded in practice, means the
same thing in every year, and is the form a conditioning feature usually wants for the
reason `08_garch_volatility` gave about ranks. The figure draws it smoothed over a
month, because the unsmoothed daily ratio is dominated by whichever single session
sits in its numerator and reaches several times the average.

```python
term_structure = (
    spy_volatility.with_columns(
        short=pl.col("volatility").rolling_mean(HORIZONS["daily"]),
        long=pl.col("volatility").rolling_mean(HORIZONS["monthly"]),
    )
    .with_columns(ratio=pl.col("short") / pl.col("long"))
    .drop_nulls()
)

frame = term_structure.select(["timestamp", "ratio"]).to_pandas().set_index("timestamp")
smoothed = frame["ratio"].rolling(HORIZONS["monthly"]).mean()

fig, ax = plt.subplots(figsize=FIGSIZE["single"])
ax.plot(frame.index, frame["ratio"], linewidth=0.4, alpha=0.4, color=COLORS["recede"])
ax.plot(frame.index, smoothed, linewidth=1, color=COLORS["blue"])
ax.axhline(1.0, color=COLORS["neutral"], linestyle="--", linewidth=0.6)
ax.fill_between(
    frame.index, 1.0, smoothed, where=smoothed > 1.0, alpha=0.25, color=COLORS["copper"]
)
ax.fill_between(frame.index, 1.0, smoothed, where=smoothed <= 1.0, alpha=0.25, color=COLORS["blue"])
ax.set_ylabel("Short horizon over long")
ax.set_xlabel("Session")
ax.set_title("Volatility above its own month, and below it")
show_with_alt(
    fig,
    "The ratio of a one-session volatility to a twenty-two-session average, drawn faintly "
    "day by day with a smoothed line over it and a dashed reference at one. The smoothed "
    "line oscillates around the reference, shaded on one side when above and the other "
    "when below, with the largest excursions above during the sharpest sell-offs.",
)
```

# Part 3: how jagged is the path

The **Hurst exponent** $H$ describes how the size of a series' movement scales with the
length of the interval measured. For a random walk the distance covered grows with the
square root of the interval, which is $H = 1/2$. Above a half the series is persistent:
a move is more likely to be followed by another in the same direction, so the path
wanders further than a random walk would. Below a half it is anti-persistent, or
**rough**: moves tend to be reversed, so the path is jagged and covers less ground.

Two things are worth separating before estimating it. Applied to **returns**, $H$ says
whether a price is trending or reverting, and for a liquid market it comes out near a
half because anything else would be a tradable pattern. Applied to **log-volatility**,
it says how jagged the volatility path is, and Gatheral, Jaisson and Rosenbaum (2018)
found it comes out near a tenth across every asset they measured, far below the half
that the standard volatility models had assumed for decades.

## Two ways to estimate it

**Rescaled range** is the original method. Cut the series into windows of size $s$;
inside each, take the cumulative deviations from that window's mean, measure the range
they cover, and divide by the window's standard deviation. Average over windows, repeat
for several sizes, and the slope of the average against the size on a log scale is $H$.

**Detrended fluctuation analysis** does the same and is unaffected by a trend. It
integrates the series first, removes a straight line from each window rather than a
mean, and measures what is left. The slope of that against the window size is the same
exponent, and it does not move when the series has a drift the rescaled range would
absorb into its estimate.

Both fit a straight line on a log scale, so both report an $R^2$ alongside the exponent,
and that $R^2$ is the check on whether the scaling relationship the method assumes holds
at all. An exponent read off a set of points that are not on a line is not an estimate
of anything.

```python
MINIMUM_WINDOW = 10
WINDOW_GROWTH = 1.5


def window_sizes(n: int, minimum: int, maximum: int | None = None) -> list[int]:
    """Geometrically spaced window sizes, so the log axis is evenly covered."""
    maximum = maximum if maximum is not None else n // 4
    sizes, size = [], minimum
    while size <= maximum:
        sizes.append(size)
        size = int(size * WINDOW_GROWTH)
    return sorted(set(sizes))


def rescaled_range(series: np.ndarray, minimum: int = 20, maximum: int | None = None):
    """Hurst exponent by rescaled range, with the points the line was fitted to."""
    log_sizes, log_statistic = [], []
    for size in window_sizes(len(series), minimum, maximum):
        values = []
        for index in range(len(series) // size):
            window = series[index * size : (index + 1) * size]
            spread = window.std(ddof=1)
            if spread > 0:
                deviations = np.cumsum(window - window.mean())
                values.append((deviations.max() - deviations.min()) / spread)
        if values:
            log_sizes.append(np.log(size))
            log_statistic.append(np.log(np.mean(values)))
    fit = stats.linregress(log_sizes, log_statistic)
    return np.array(log_sizes), np.array(log_statistic), fit.slope, fit.rvalue**2


def detrended_fluctuation(
    series: np.ndarray, minimum: int = MINIMUM_WINDOW, maximum: int | None = None, order: int = 1
):
    """Hurst exponent by detrended fluctuation analysis, with its fitted points."""
    profile = np.cumsum(series - series.mean())
    log_sizes, log_fluctuation = [], []
    for size in window_sizes(len(series), minimum, maximum):
        steps = np.arange(size)
        fluctuations = []
        for index in range(len(series) // size):
            segment = profile[index * size : (index + 1) * size]
            trend = np.polyval(np.polyfit(steps, segment, order), steps)
            fluctuations.append(np.sqrt(np.mean((segment - trend) ** 2)))
        if fluctuations:
            log_sizes.append(np.log(size))
            log_fluctuation.append(np.log(np.mean(fluctuations)))
    fit = stats.linregress(log_sizes, log_fluctuation)
    return np.array(log_sizes), np.array(log_fluctuation), fit.slope, fit.rvalue**2
```

## The two series, and what each one says

Returns go in directly. Volatility goes in as the **increments of its logarithm**,
which is the quantity the rough-volatility result is about: the model it argues against
says those increments behave like a random walk's, and the finding is that they do not.

```python
returns_series = spy_volatility["log_return"].to_numpy()
log_volatility = np.log(spy_volatility["volatility"].to_numpy().clip(min=1e-10))
log_volatility_increments = np.diff(log_volatility)

estimates = []
for name, series in [
    ("returns", returns_series),
    ("log-volatility increments", log_volatility_increments),
]:
    rs_sizes, rs_values, rs_exponent, rs_fit = rescaled_range(series)
    dfa_sizes, dfa_values, dfa_exponent, dfa_fit = detrended_fluctuation(series)
    estimates.append(
        {
            "series": name,
            "H, rescaled range": rs_exponent,
            "R-squared, rescaled range": rs_fit,
            "H, detrended fluctuation": dfa_exponent,
            "R-squared, detrended fluctuation": dfa_fit,
            "_points": (rs_sizes, rs_values, dfa_sizes, dfa_values),
        }
    )

display(pd.DataFrame(estimates).drop(columns="_points"))
```

```python
fig, axes = plt.subplots(2, 2, figsize=FIGSIZE["grid_3x2"])

for row, estimate in enumerate(estimates):
    rs_sizes, rs_values, dfa_sizes, dfa_values = estimate["_points"]
    for column, (sizes, values, exponent, label) in enumerate(
        [
            (rs_sizes, rs_values, estimate["H, rescaled range"], "Rescaled range"),
            (dfa_sizes, dfa_values, estimate["H, detrended fluctuation"], "Detrended fluctuation"),
        ]
    ):
        ax = axes[row, column]
        ax.scatter(sizes, values, s=14, color=COLORS["blue"], zorder=3)
        ax.plot(
            sizes,
            exponent * sizes + (values[0] - exponent * sizes[0]),
            linestyle="--",
            linewidth=1,
            color=COLORS["negative"],
        )
        ax.set_title(f"{label}, {estimate['series']}")
        ax.set_xlabel("Log window size")
        ax.set_ylabel("Log statistic")

fig.suptitle("The exponent is a slope, so the points have to be on a line")
show_with_alt(
    fig,
    "A two by two grid of log-log scatter plots with a dashed fitted line through each. "
    "The top row is for returns and the bottom for log-volatility increments; the left "
    "column is the rescaled range and the right the detrended fluctuation. The returns "
    "points lie close to their lines with a clear positive slope; the log-volatility "
    "points lie on much flatter lines.",
)
```

```python
for estimate in estimates:
    print(
        f"{estimate['series']}: "
        f"rescaled range {estimate['H, rescaled range']:.3f}, "
        f"detrended fluctuation {estimate['H, detrended fluctuation']:.3f}, "
        f"against 0.5 for a random walk"
    )
```

Read the two rows against a half and against each other. Whichever method is used, the
two series land in different places, and the direction of the difference is the finding:
returns scale close to a random walk's, and the increments of log-volatility scale far
below it, meaning volatility reverses direction much more often than a random walk does.

That is what "rough" means, and its consequence is practical. A model built on the
assumption that log-volatility is a random walk, which describes most of the standard
volatility models, is assuming a path smoother than the one the data traces. It will
be slow to follow volatility down after a spike, because it expects the level to persist
when in fact the next increment is more likely than not to reverse the last.

The two methods do not agree exactly, and neither is a measurement of a physical
constant. Each estimates a slope from a handful of points, the points are averages over
overlapping structure, and the estimate moves with the window range chosen. Read the
exponent together with the $R^2$ beside it, and treat a difference of a few hundredths
between methods as noise rather than as a finding.

## Roughness over time

Estimated on a moving window, the exponent becomes a column. For returns it is a slow
statement about whether the market is trending or reverting; for log-volatility it says
how bursty the current volatility regime is.

```python
HURST_WINDOW = 252
HURST_STEP = 10  # sessions between estimates; each estimate is a fit over many sub-windows

rolling_rows = []
for end in range(HURST_WINDOW, len(returns_series), HURST_STEP):
    _, _, returns_exponent, _ = detrended_fluctuation(
        returns_series[end - HURST_WINDOW : end], maximum=HURST_WINDOW // 4
    )
    # The increments series is one shorter than the returns series, so its window ends one
    # index earlier and covers the same sessions.
    _, _, volatility_exponent, _ = detrended_fluctuation(
        log_volatility_increments[max(0, end - HURST_WINDOW - 1) : end - 1],
        maximum=HURST_WINDOW // 4,
    )
    rolling_rows.append(
        {
            "timestamp": spy_volatility["timestamp"][end],
            "returns": returns_exponent,
            "log-volatility": volatility_exponent,
        }
    )

rolling_hurst = pd.DataFrame(rolling_rows).set_index("timestamp")
print(f"Rolling estimates: {len(rolling_hurst)}, each fitted on {HURST_WINDOW} sessions")
print(
    "Median exponent: "
    f"returns {rolling_hurst['returns'].median():.3f}, "
    f"log-volatility {rolling_hurst['log-volatility'].median():.3f}"
)
print(
    "Share of windows above one half: "
    f"returns {(rolling_hurst['returns'] > 0.5).mean():.1%}, "
    f"log-volatility {(rolling_hurst['log-volatility'] > 0.5).mean():.1%}"
)
```

```python
fig, axes = plt.subplots(2, 1, figsize=FIGSIZE["dual_v"], sharex=True)

for ax, column, color, title in [
    (axes[0], "returns", COLORS["blue"], "Returns: trending above the line, reverting below"),
    (axes[1], "log-volatility", COLORS["copper"], "Log-volatility: rough throughout"),
]:
    ax.plot(rolling_hurst.index, rolling_hurst[column], linewidth=0.9, color=color)
    ax.axhline(0.5, color=COLORS["negative"], linestyle="--", linewidth=0.7)
    ax.annotate(
        "A random walk",
        xy=(0.995, 0.5),
        xycoords=("axes fraction", "data"),
        xytext=(0, 3),
        textcoords="offset points",
        ha="right",
        fontsize=7,
        color=COLORS["negative"],
    )
    ax.set_ylabel("Exponent")
    ax.set_title(title)
axes[1].set_xlabel("Session the window ends")

fig.suptitle(f"The exponent on a {HURST_WINDOW}-session moving window")
show_with_alt(
    fig,
    "Two stacked panels, each drawing a rolling Hurst exponent against a dashed reference "
    "at one half. The top, for returns, oscillates around the reference and crosses it "
    "repeatedly. The bottom, for log-volatility increments, stays well below the reference "
    "throughout the sample without approaching it.",
)
```

The top panel is a feature and the bottom panel is a fact about markets. Returns cross
the line in both directions, so a column holding that exponent separates periods; the
log-volatility exponent never comes near it, so a column holding that one separates
almost nothing and its value is a property of the asset class rather than of the moment.

That asymmetry is the practical reading of the rough-volatility literature. The finding
is important because it says the standard model is wrong, and it is a poor feature for
exactly the same reason: it is wrong everywhere, in the same direction, all the time.

## The library expression, and what it measures

`hurst_exponent` computes a rolling estimate as a Polars expression, and it takes a
**price** column. That is a third series, different from both of the ones above: the
rescaled range of a price level rather than of its returns or of its volatility. The
three answer different questions and there is no reason for them to agree.

```python
library_hurst = spy.with_columns(hurst=hurst_exponent("close", period=HURST_WINDOW))
library_values = library_hurst["hurst"].drop_nulls()

display(
    pd.DataFrame(
        [
            {
                "series": "close price, library expression",
                "median": library_values.median(),
                "share above one half": float((library_values > 0.5).mean()),
            },
            {
                "series": "returns, detrended fluctuation",
                "median": rolling_hurst["returns"].median(),
                "share above one half": float((rolling_hurst["returns"] > 0.5).mean()),
            },
        ]
    )
)
```

The two differ, and the difference is the input rather than the implementation. Use the
expression where a rolling exponent of the price is what is wanted and the cost of the
manual version across a panel is prohibitive; use the manual version where the series
being characterized is a return or a volatility, which is most of the time in this
chapter.

# Part 4: across the panel

Every number above comes from one symbol. Fitting the same three quantities across the
ETF panel says which of them are properties of this symbol and which are properties of
the asset class.

```python
MINIMUM_PANEL_SESSIONS = 600

panel_source = (
    load_etfs()
    .filter(pl.col("timestamp") >= pl.lit(START_DATE).str.to_date())
    .filter(pl.col("timestamp") <= pl.lit(END_DATE).str.to_date())
    .sort(["symbol", "timestamp"])
)


def panel_row(symbol: str) -> dict | None:
    """HAR coefficients and both exponents for one symbol, or None where the history is short."""
    frame = panel_source.filter(pl.col("symbol") == symbol).select(
        ["timestamp", "open", "high", "low", "close"]
    )
    if frame.height < MINIMUM_PANEL_SESSIONS:
        return None

    volatility = range_volatility(frame)
    model_frame = har_frame(volatility)
    if model_frame.height < MINIMUM_PANEL_SESSIONS // 2:
        return None

    fit = fit_har(model_frame)
    increments = np.diff(np.log(volatility["volatility"].to_numpy().clip(min=1e-10)))
    _, _, returns_exponent, _ = detrended_fluctuation(volatility["log_return"].to_numpy())
    _, _, volatility_exponent, _ = detrended_fluctuation(increments)

    return {
        "symbol": symbol,
        **{name: fit.params[name] for name in HORIZONS},
        "r_squared": fit.rsquared,
        "H returns": returns_exponent,
        "H log-volatility": volatility_exponent,
    }


panel = pd.DataFrame(
    [row for symbol in panel_source["symbol"].unique().sort() if (row := panel_row(symbol))]
)
print(f"Symbols fitted: {len(panel)} of {panel_source['symbol'].n_unique()}")
display(panel.describe().loc[["mean", "std", "min", "max"]].round(4))
```

```python
fig, axes = plt.subplots(2, 2, figsize=FIGSIZE["grid_3x2"])

ax = axes[0, 0]
for name, color in zip(HORIZONS, [COLORS["blue"], COLORS["amber"], COLORS["copper"]]):
    ax.hist(panel[name], bins=25, alpha=0.6, color=color, label=name)
ax.axvline(0, color=COLORS["neutral"], linestyle="--", linewidth=0.6)
ax.set_xlabel("Coefficient")
ax.set_ylabel("Symbols")
ax.set_title("The three HAR weights across the panel")
ax.legend(fontsize=7)

ax = axes[0, 1]
ax.hist(panel["r_squared"], bins=25, color=COLORS["blue"])
ax.set_xlabel("R-squared")
ax.set_ylabel("Symbols")
ax.set_title("How much of volatility the model explains")

for ax, column, reference, title in [
    (axes[1, 0], "H returns", 0.5, "Exponent on returns"),
    (axes[1, 1], "H log-volatility", 0.5, "Exponent on log-volatility"),
]:
    ax.hist(panel[column], bins=25, color=COLORS["copper"])
    ax.axvline(reference, color=COLORS["negative"], linestyle="--", linewidth=0.8)
    ax.set_xlabel("Exponent")
    ax.set_ylabel("Symbols")
    ax.set_title(title)

fig.suptitle("What varies across symbols, and what does not")
show_with_alt(
    fig,
    "A two by two grid of histograms across the ETF panel. Top left overlays the three HAR "
    "coefficients, which occupy different and partly overlapping ranges. Top right is the "
    "R-squared, concentrated over a moderate range. Bottom left is the exponent on returns, "
    "clustered near the dashed reference at one half. Bottom right is the exponent on "
    "log-volatility, clustered far below that reference with no overlap between the two.",
)
```

```python
print(
    "Exponent on returns: "
    f"median {panel['H returns'].median():.3f}, "
    f"{(panel['H returns'] > 0.5).sum()} of {len(panel)} symbols above one half"
)
print(
    "Exponent on log-volatility: "
    f"median {panel['H log-volatility'].median():.3f}, "
    f"{(panel['H log-volatility'] > 0.5).sum()} of {len(panel)} symbols above one half"
)
print(
    "HAR coefficient with the largest median: "
    + max(HORIZONS, key=lambda name: panel[name].median())
)
```

The two bottom panels are the answer to the question the panel was fitted to settle. The
exponent on returns clusters close to a half with symbols on both sides of it, which is
what a liquid market should look like. The exponent on log-volatility sits far below a
half for every symbol, with no overlap at all between the two distributions, which
reproduces the rough-volatility finding across a hundred assets rather than one.

Compare the returns figure here against the rolling one two sections up before reading
either as a level. This one estimates the exponent once over each symbol's whole history;
that one estimates it over 252 sessions at a time, and their medians differ. An exponent
is a slope through averaged points and the estimate moves with the range of window sizes
the fit was given, so the two numbers are answers to slightly different questions rather
than a disagreement about the market. What both agree on is the direction, which is that
returns sit near a half and volatility does not.

The HAR coefficients vary more across symbols than either exponent does, and that
variation is what makes them worth carrying per symbol rather than fitting once.

## The features this notebook produces

| Column | What it is | Causal |
|---|---|---|
| `daily`, `weekly`, `monthly` | lagged averages of range-based volatility over three horizons | yes |
| HAR coefficients | the weight on each horizon, refitted on a moving window | yes, refit on past windows |
| term structure ratio | the short horizon over the long one | yes |
| exponent on returns | how trending or reverting the price has been, on a moving window | yes |
| exponent on log-volatility | how jagged the volatility path has been | yes |

The realized volatility measured from minute bars is not on the list, because it exists
for two symbols and two years and the point of Part 1 was to establish which daily
estimator stands in for it.

## Key takeaways

1. **Check the estimator against a measurement where you can.** A squared close-to-close
   return misses a session that moved and came back; the range-based estimators see it,
   and comparing them against volatility measured from minute bars is how you find out by
   how much and in which direction each is biased.
2. **HAR is a linear regression, and that is its advantage.** Three lagged averages track
   realized volatility about as closely as a fitted GARCH does, and unlike GARCH each
   coefficient is readable as the weight one horizon carries.
3. **Overlapping regressors are a reason to widen the standard errors, not a proof the
   residuals are correlated.** A correctly specified model can have independent
   innovations despite rolling-average regressors. To read what one relaxed assumption
   costs, change only that assumption: the same Newey-West estimator at zero lags and at
   the model's longest horizon differs in the window alone, which no comparison across two
   different covariance estimators can claim.
4. **Two models forecasting differently-measured targets cannot be ranked by their
   errors.** Report the mean of each forecast next to the mean of the target, and read
   the correlation, which the level difference does not touch.
5. **Roughness is a finding about markets and a poor feature.** The exponent on
   log-volatility sits far below a half for every symbol in the panel, which is why the
   standard models are wrong and also why a column holding it separates nothing. The
   exponent on returns crosses a half in both directions and is the one worth carrying.

**Known limitations.** Part 1 compares estimators on one symbol over two years, both
chosen by what intraday data exists rather than by what would be representative. The HAR
and GARCH comparison uses one split at one date. Every Hurst estimate is a slope through
a handful of averaged points and moves with the window range chosen, so differences of a
few hundredths between methods carry no information. And the panel's ETFs overlap
heavily in what they hold, so its hundred rows are far fewer than a hundred independent
observations.

**Previous**: `08_garch_volatility` for a single-horizon volatility model.
**Next**: `10_uncertainty_features`, which treats the model's own uncertainty as the
feature rather than its estimate.
![notebook output](figures/p1_1.png)
![notebook output](figures/p1_2.png)
![notebook output](figures/p1_3.png)
![notebook output](figures/p1_4.png)
![notebook output](figures/p1_5.png)
![notebook output](figures/p1_6.png)
![notebook output](figures/p1_7.png)
![notebook output](figures/p1_8.png)