Comparing Regression Losses for Short-Straddle Return Ranking

# Short straddles: three loss functions against one rank metric


# Short straddles: three loss functions against one rank metric

[`06_linear`](06_linear.ipynb) set out the problem this case study exists to teach. A short
at-the-money straddle can earn at most the premium it collected and can lose without bound, so
its return distribution is capped on one side and open on the other. A linear regression fits by
minimizing squared error, which weights an observation by the square of its distance from the
prediction; the information coefficient is a rank correlation, which treats a loss of twenty
percent and a loss of nine thousand percent as adjacent. On this target those two disagree, and
the linear notebook could only pose that as a question, because ordinary least squares has one
loss function and no alternative to compare it against.

**This notebook answers it.** The declared grid fits every capacity setting under three
objectives - squared error, absolute error, and Huber, which behaves like squared error near the
middle of the distribution and like absolute error out in the tail. Nothing else changes: the
same features, the same folds, the same rows, the same metric. If the mismatch between the
fitting loss and the reporting metric is what holds this family back, the objectives should
separate, and they should separate in the order of how much weight each one gives the extremes.

Two further dials come with the model class:

- **Capacity**, set by `num_leaves`: how finely one tree may partition the feature space, and so
  how sharply it can separate the conditions under which the position keeps its premium from
  those where it does not.
- **When to stop**, set by the number of trees. A boosted model has a meaningful state at every
  iteration, so each configuration is scored at ten points along its own training run. Fifteen
  declared configurations at ten checkpoints is 150 candidate models, and **a checkpoint is part
  of a configuration rather than a detail of how it was fitted.**

Two folds, validating on 2019 and 2020 - the second a year in which a short-volatility position
experienced the event it exists to be paid for. Every number below is an average over those two
periods, and the second dominates the tail behaviour the objectives are being compared on.

**Learning objectives.** By the end of this notebook you will be able to:

- Say what a loss function claims about which errors matter, and check that claim against the
  metric a result is judged on.
- Read a comparison of squared error, absolute error and Huber loss on a target with a hard
  ceiling and an unbounded downside, and say what their ordering means.
- Explain why a boosted model produces one result per checkpoint while a linear model produces
  one result in total, and what that implies for counting candidates.
- Read a learning curve of out-of-sample information coefficient against tree count, and say
  whether an apparent peak is a turning point or the highest of ten noisy readings.
- Judge a result by the number of candidates behind it and the amount of validation data under
  it, rather than by its rank in a table.

**Book reference**: Chapter 12, Section 12.2 (GBM libraries) and Section 12.3 (how to tune a
boosted model). Chapter 6, Section 6.7 (Search accounting and run logging) introduces the run
log this notebook writes to.

**Prerequisites**: [`03_financial_features`](03_financial_features.ipynb) and
[`04_model_based_features`](04_model_based_features.ipynb) have written the feature matrices,
[`05_evaluation`](05_evaluation.ipynb) has established the walk-forward folds, and
[`06_linear`](06_linear.ipynb) fitted the linear population this one is compared against.

**What it writes**: one training run per configuration and one complete validation prediction
set per configuration and checkpoint, in `run_log/registry.db` and under `run_log/training/` and
`run_log/predictions/`, grouped under a named population.
[`12_backtest`](12_backtest.ipynb) reads that population and selects on validation backtest
Sharpe. **Selection happens there, not here.**

```python
"""Fit the declared short-straddle gradient boosting population on the walk-forward folds."""

import plotly.graph_objects as go
import polars as pl

from case_studies.research import (
    declared_labels,
    load_model_configs,
    model_requests,
    narrows_declared_catalog,
    open_study,
    resolved_model_plan,
    run_model_population,
    supersedes_for_run,
)
from utils.style import COLORS, show_plotly_with_alt
```

```python
LABELS: list[str] = []
EXECUTION_TIER = "canonical"
WORKSPACE: str = ""
PREVIEW_REDUCTIONS: dict = {}
CONFIG_NAMES: list[str] = []
POPULATION_NAME = ""
SUPERSEDES_POPULATION: str = ""
```

```python
study = open_study("sp500_options", execution_tier=EXECUTION_TIER, workspace=WORKSPACE or None)
```

## 1. Which labels, and which models

The labels are the same set the linear notebook fitted - the ones `config/setup.yaml` puts in
the sweep, which here is `ret_to_expiry` alone: the return the short straddle earns from the
decision date to expiry. Keeping the set fixed is what makes the two populations comparable;
the families differ, the targets do not. `06_linear` says why `config/training/` holds five
menus when the sweep declares one label.

```python
declared_labels(study, "gbm")
```

Each declared label's menu at `config/training/{label}.yaml` lists 15 named configurations
under `gbm:`, and each resolves to a preset in `case_studies/config/lgb/`. The grid is a
product of two axes:

- **Five capacity profiles.** `default` uses the library's own leaf count; the rest fix it at 7,
  15, 31 and 63. Leaf count is the direct control on how finely one tree may partition the
  feature space. It is the dial that decides how finely a model can separate the conditions
  under which a short straddle keeps its premium from the ones where it does not.
- **Three objectives.** `mse` minimizes squared error, `mae` absolute error, and `huber` behaves
  like squared error for small residuals and like absolute error beyond a threshold derived from
  each fold's own label spread.

Every configuration runs the same number of boosting iterations with the same learning rate, so
the grid isolates capacity and loss rather than confounding them with training length.

```python
configs = load_model_configs(
    study,
    "gbm",
    labels=LABELS or None,
    config_names=CONFIG_NAMES or None,
)
# `model_class` is empty for every row: the LightGBM presets under `case_studies/config/lgb/`
# declare their objective in `params` and carry no `model_class` key, unlike the linear presets
# where it names the estimator. Showing the column would put a blank field in front of the reader.
configs.drop("model_class")
```

`LABELS` and `CONFIG_NAMES` both narrow what is fitted, and a narrowed run declares a different
set of members than the canonical population does. A population is immutable once written, so
such a run must publish under its own name, and this says so here rather than several cells
later in a message about hashes.

```python
if narrows_declared_catalog(study, "gbm", configs) and not POPULATION_NAME:
    raise ValueError(
        f"this run declares {configs.height} label-configuration pairs, which is not the "
        f"complete declared catalog, so it cannot publish the canonical population; pass "
        f"POPULATION_NAME to give it its own"
    )
```

## 2. Binding the declarations to the data

Resolving reads the label and feature files, computes the fold boundaries, works out the exact
rows each fit must predict, and turns any data-dependent parameter into the number it will use.
Huber's threshold is one of those: it is a fraction of the training labels' standard deviation,
so it is a different number on every fold and is resolved from that fold's own data.

Nothing is fitted here, so the plan can be inspected first. Four things to check:

- **`feature_count`, `eligible_entities` and `eligible_rows` agree across every row.** A row that
  differs is a configuration measured on a different sample from its neighbours.
- **`folds` is the same everywhere**, and equals the number of walk-forward splits.
- **`validation_start` and `validation_end` bracket the development sample**, with none of the
  held-out tail visible.
- **`checkpoints` is where this differs from the linear plan.** It is the number of training
  states each configuration will publish predictions for. Multiply it by the number of rows to
  get the number of candidate models this notebook is about to create.

```python
requests = model_requests(
    study,
    configs,
    execution_tier=EXECUTION_TIER,
    preview_reductions=PREVIEW_REDUCTIONS,
)
resolved = tuple(request.resolve() for request in requests)

plan = resolved_model_plan(resolved)
plan.select(
    "config_name",
    "feature_count",
    "eligible_entities",
    "eligible_rows",
    "folds",
    "checkpoints",
    "validation_start",
    "validation_end",
)
```

## 3. Fitting the population

`run_model_population` fits every resolved request. For one request it walks the folds, and on
each one:

1. takes the rows inside that fold's training window,
2. casts the design matrix to the precision LightGBM works in and leaves missing values in
   place - a tree routes a missing value down its own branch, so imputing a median here would
   hand the model an observation nobody made,
3. fits the declared number of boosting iterations,
4. predicts the fold's validation rows at each checkpoint, using only the trees built up to that
   iteration.

Step 4 is what makes one fit produce many results. The fold predictions are concatenated into
one series per checkpoint covering the whole validation period, and each becomes its own
registered prediction set with its own identity.

There are two ways to submit these fits and this notebook takes the first of them. Resolving
every request before submitting is what let the plan table above show real feature counts, fold
geometry and eligible-row counts, and it costs one fold preparation per configuration: each
resolved request holds its own prepared folds and runs on its own. Submitting the requests
unresolved instead reaches a family batch runner that walks folds on the outside and
configurations on the inside, so one prepared fold set serves every configuration - at the cost
of a plan built from placeholder folds rather than the real ones. The choice is between seeing
what will be fitted and holding less at once; on a grid this size the first is affordable.

**What the call publishes is a population**: a named, immutable list of the prediction sets it
will produce, written down before the first fit. Afterwards every member must exist and be
complete, which is what makes the downstream comparison well defined.

The default name is the contract with `11_model_analysis` and `12_backtest`, which resolve this
population by name, rather than a label of convenience; a run that narrows the member set has to
pass its own.

`SUPERSEDES_POPULATION` names the earlier snapshot this one replaces. A population is hashed
over its members *and* over what it supersedes, so the default carries the hash the published
snapshot actually superseded: leaving it empty would compute a different population and the
registry would refuse it against the one on record.

The default hash records what the published snapshot replaced, so it is only correct under that
population's name and in a registry that holds that generation. `supersedes_for_run` is that
decision: it withholds the hash for a preview, whose population lives and dies with its
workspace and has nothing to extend, for a run under a name of its own, whose first version
would be refused for claiming a predecessor, and for a reader's clean clone, which holds no
generation at all. The comparison this replaced tested the declaration against a second
hardcoded hash instead of against the registry, so a custom name kept a predecessor that name
never had.

```python
DEFAULT_POPULATION_NAME = "sp500-options-gbm-validation-v1"

population_name = POPULATION_NAME or DEFAULT_POPULATION_NAME
supersedes = supersedes_for_run(
    study,
    population_name=population_name,
    declared=SUPERSEDES_POPULATION or None,
    execution_tier=EXECUTION_TIER,
)
execution, population = run_model_population(
    study,
    resolved,
    population_name=population_name,
    supersedes=supersedes,
)

fitted = sum(len(item["fitted_folds"]) for item in execution.diagnostics)
reused = sum(len(item["reused_folds"]) for item in execution.diagnostics)
print(f"{len(execution.runs)} configurations: {fitted} folds fitted, {reused} reused")
print(f"population {population.name}: {len(population.members)} prediction sets")
```

Re-running this notebook unchanged costs the time it takes to read the data. Every identity is
re-derived from the inputs, the registry already holds the matching rows, and the runner returns
the stored result rather than fitting again.

### Running configurations of your own

The published run log is read-only. To add runs, open the study against a workspace, which holds
its own registry and artifacts and reads the same labels and features:

```python
study = open_study("sp500_options", workspace="~/ml4t-experiments")
configs = load_model_configs(
    study, "gbm", labels=["ret_to_expiry"], config_names=["leaves_15_huber", "leaves_31_huber"]
)
requests = model_requests(study, configs)
resolved = tuple(request.resolve() for request in requests)
execution, population = run_model_population(study, resolved, population_name="my-gbm-v1")
```

`CONFIG_NAMES` fits a subset of what the menu declares. To fit something new, add a preset at
`case_studies/config/lgb/leaves_127_huber.yaml` and list `leaves_127_huber` under `gbm:` in the
label's menu. Editing an existing preset changes that configuration's identity, so its result
registers as a new row beside the old one rather than replacing it.
[`RUN_LOG.md`](../RUN_LOG.md#running-your-own-configurations) covers the rest.

## 4. What came out

One row per configuration and checkpoint. `ic_mean` is the **information coefficient**: on each
validation date, rank the positions by the model's prediction, rank them by the return they went
on to earn, correlate the two rankings, and average that daily correlation over the validation
period.

The table is sorted by IC, and the top of it is the trap this notebook exists to describe. The
leading row is the maximum of 150 numbers. Reading it as the result of one experiment would
attribute to the model whatever the stopping point contributed, and the section below measures
how large that contribution is before anything is concluded from the ranking.

```python
catalog = execution.catalog_rows.select(
    "config_name",
    "label",
    "complete",
    "checkpoint_value",
    "ic_mean",
    "ic_std",
    "ic_n_days",
    "n_folds",
    "training_hash",
    "prediction_hash",
).sort("ic_mean", descending=True)

if catalog.filter(~pl.col("complete")).height:
    raise RuntimeError("gbm execution returned a partial prediction set")
```

A configuration scores on every date its folds cover, and a configuration that covers fewer
dates than its neighbours is measured over a shorter window, so its IC is not comparable with
theirs. Coverage is judged against each label's own maximum rather than a global one: with one
label in the sweep the two are the same number, and adding a variant to `setup.yaml` is all it
would take for a global maximum to mark a whole grid incomplete for a reason that has nothing
to do with the models.

The charts below are one panel, which is right only while the sweep is one label, so a second
label stops the notebook rather than pooling two horizons into one ranking.

```python
catalog = catalog.with_columns(
    full_coverage=pl.col("ic_n_days") == pl.col("ic_n_days").max().over("label")
)
if catalog.get_column("label").n_unique() > 1:
    raise NotImplementedError(
        "this notebook charts one label; facet the figures before adding a sweep variant"
    )

print(f"{catalog.height} candidate models: {catalog.n_unique('config_name')} configurations")
print(f"at {catalog.n_unique('checkpoint_value')} checkpoints each")
catalog.select(
    "config_name",
    "checkpoint_value",
    "ic_mean",
    "ic_std",
    "ic_n_days",
    "full_coverage",
).head(15)
```

### What more trees do

Each line traces one configuration's out-of-sample IC as trees are added to it. This is the
figure the checkpoint dimension exists to produce, and it separates two things a single
end-of-training number cannot.

A line that rises and then falls has an interior optimum: the model was still learning, then
began fitting the training window at the expense of the validation folds. A line that wanders
without trend around zero never had anything to learn in the first place, and its highest point
is wherever the noise happened to peak. The difference matters, because both produce a
respectable-looking maximum.

```python
curves = catalog.filter("full_coverage").sort("config_name", "checkpoint_value")
objectives = {"mse": COLORS["blue"], "mae": COLORS["amber"], "huber": COLORS["copper"]}


def objective_of(name: str) -> str:
    """Read the loss function out of a declared configuration name."""
    return next((key for key in objectives if name.endswith(key)), "mse")
```

One trace per configuration, coloured by the objective its name declares, so the question of
whether the loss function is what separates them can be read off the chart rather than argued.

```python
fig_curves = go.Figure()
for config_name in curves.get_column("config_name").unique(maintain_order=True):
    series = curves.filter(pl.col("config_name") == config_name)
    fig_curves.add_trace(
        go.Scatter(
            x=series.get_column("checkpoint_value").to_list(),
            y=series.get_column("ic_mean").to_list(),
            mode="lines",
            name=config_name,
            line=dict(color=objectives[objective_of(config_name)], width=1.5),
            opacity=0.75,
        )
    )
fig_curves.add_hline(y=0, line_width=1, line_dash="dash", line_color=COLORS["neutral"])
fig_curves.update_layout(
    title="Validation IC against boosting iteration, by loss function",
    height=550,
    width=1000,
    margin=dict(t=70),
    legend=dict(font=dict(size=9)),
)
fig_curves.update_xaxes(title_text="Boosting iterations (trees kept)")
fig_curves.update_yaxes(title_text="Mean cross-sectional IC (validation)")
show_plotly_with_alt(
    fig_curves,
    "Line chart of mean validation information coefficient against boosting iteration, one line "
    "per configuration, coloured by loss function: blue for squared error, amber for absolute "
    "error, copper for Huber. A dashed line marks zero. The blue squared-error lines run below "
    "it across the whole range without exception. The copper Huber and amber absolute-error "
    "lines run above it for most of the range and are interleaved with one another rather than "
    "separated; two Huber lines and one absolute-error line dip below zero as trees are added. "
    "Ten of the fifteen reach their highest point at the first or second checkpoint and drift "
    "afterwards.",
)
```

### Whether the loss function is what separates them

The curves are coloured by objective because that is the axis with a mechanism behind it. If
heavy tails are steering the squared-error fits, the three colours should separate, and they
should separate more as trees are added, since each additional tree is fitted to the residuals
the previous ones left.

The chart below drops the checkpoint dimension by taking each configuration's final state, so
every configuration is compared at the same amount of training. That is the comparison that does
not require choosing anything after the fact.

```python
final_iteration = int(catalog.get_column("checkpoint_value").max())
final = (
    catalog.filter(pl.col("checkpoint_value") == final_iteration)
    .filter("full_coverage")
    .with_columns(objective=pl.col("config_name").map_elements(objective_of, return_dtype=pl.Utf8))
    .sort("ic_mean", descending=True)
)

fig_obj = go.Figure(
    go.Bar(
        x=final.get_column("config_name").to_list(),
        y=final.get_column("ic_mean").to_list(),
        marker_color=[objectives[value] for value in final.get_column("objective")],
        text=[f"{value:+.3f}" for value in final.get_column("ic_mean")],
        textposition="outside",
        cliponaxis=False,
    )
)
fig_obj.add_hline(y=0, line_width=1, line_dash="dash", line_color=COLORS["neutral"])
fig_obj.update_layout(
    title="Validation IC at the final iteration, coloured by loss function",
    height=500,
    width=1000,
    showlegend=False,
    margin=dict(t=70),
)
fig_obj.update_xaxes(title_text="Configuration (sorted by validation IC)", tickangle=-45)
fig_obj.update_yaxes(title_text="Mean cross-sectional IC (validation)")
show_plotly_with_alt(
    fig_obj,
    "Bar chart of mean validation information coefficient for every full-coverage configuration "
    "at its final boosting iteration, sorted descending and coloured by loss function, against a "
    "dashed zero line. The bars divide by colour rather than by leaf count: every blue "
    "squared-error bar hangs below the line, and together they take the five rightmost places "
    "with nothing else among them. The copper Huber and amber absolute-error bars are "
    "interleaved to the left of that block, mostly above the line, with two copper and one "
    "amber bar just below it. Leaf count varies within each colour group without ordering it.",
)
```

### How much the checkpoint moves a configuration

One number per configuration: the range its IC covers across its own ten checkpoints. This is
the quantity that decides whether choosing a stopping point is a decision worth making carefully
or one being made by noise. A configuration whose IC varies more across its own training run
than the configurations vary among themselves is one where the checkpoint, not the model, is
doing the ranking.

```python
spread = (
    curves.group_by("config_name")
    .agg(
        ic_min=pl.col("ic_mean").min(),
        ic_max=pl.col("ic_mean").max(),
        ic_final=pl.col("ic_mean").filter(pl.col("checkpoint_value") == final_iteration).first(),
    )
    .with_columns(checkpoint_range=pl.col("ic_max") - pl.col("ic_min"))
    .sort("checkpoint_range", descending=True)
)
across_configs = float(final.get_column("ic_mean").max() - final.get_column("ic_mean").min())
print(f"IC range across configurations at the final iteration: {across_configs:.4f}")
print(
    f"median IC range within one configuration: {spread.get_column('checkpoint_range').median():.4f}"
)
spread
```

## 5. What to notice

**Squared error separates from the two robust objectives, and that is the whole of the
separation.** Every squared-error configuration ends training below zero, and every one of them
ends below every Huber and every absolute-error configuration: the five squared-error settings
take the bottom five places in the ranking with nothing else among them. Above that block the
other two objectives interleave rather than order. The leading Huber and the leading absolute-
error configuration tie for first, two Huber settings end just below zero, and four of the five
absolute-error settings end above it. So what the grid supports is squared error against the
rest, not the three-way ordering the mechanism might lead you to expect. Capacity orders nothing
either way: the shallowest setting appears both in the top three and in last place.

**The mechanism is the payoff shape, and the linear result is consistent with it.** A short
straddle earns at most the premium and can lose many multiples of it, so the label's mass sits
in a narrow band with a long tail on one side. Squared error weights an observation by the square
of its error, so a handful of expiries where the position lost heavily dominate what every
successive tree is fitted to - and getting those few right does nothing for a rank correlation,
which only asks whether the positions were ordered correctly. Huber caps how much any one
observation can pull, spends its capacity on the body of the distribution where most of the
ranking is decided, and ends above squared error everywhere in this grid. The linear model in
`06_linear` minimizes squared
error and finished close to where the squared-error configurations here finish, which is what
this reading predicts.

**The lesson generalizes past this case study.** Nothing above is specific to options. It applies
whenever the target has a long tail and the metric is rank-based, which describes most of the
labels in this book. **Choosing a loss function is choosing which errors the model is paid to
reduce, and if that choice does not match the metric the result is judged on, the model will
spend its capacity where the metric gives nothing back.** The usual remedies - winsorizing the
label, ranking it, or changing the objective - are the same intervention seen from three sides.

**Read the magnitudes with the sample in mind.** The positive numbers are small, they rest on
two validation folds covering 2019 and 2020, and the second of those is one event that a
short-volatility position exists to be paid for. A single regime can move a two-fold average a
long way. What the results support is squared error's separation from the other two, which is
consistent across the grid and has a mechanism behind it; what they do not support is a ranking
of Huber against absolute error, or any particular configuration's value.

**Most configurations peak early.** Ten of the fifteen reach their highest IC at the first or
second checkpoint - inside the first fifth of the declared iterations - and drift afterwards. Combined with the objective finding, the
picture is of a model that extracts what is available quickly and then spends the remaining
iterations fitting the tail of the training window - which is precisely what squared error
rewards it for doing, and part of why that objective ends up furthest below zero.

**The checkpoint is a smaller dial than the objective here.** Across the fifteen configurations
at fixed training length the IC spans about three times what a single configuration spans across
its own ten checkpoints. That is the unusual case: in this case study the thing that decides the
answer is a declared modelling choice rather than where training happened to stop. It does not
make the leading row safe to quote - it is still the maximum of 150 numbers - but it does mean
the pattern in the second chart is the result rather than an artifact of stopping points.

**None of this selects anything.** IC measures whether predictions rank positions correctly, not
whether a strategy trading them makes money after costs. That gap is unusually wide for short
volatility: a strategy can rank correctly and still be ruined by the tail it ranked, because the
metric treats a large loss and a small one as adjacent while the profit and loss does not.
Selection is on validation backtest Sharpe over the population just published, in
[`12_backtest`](12_backtest.ipynb), where the checkpoint is part of what is selected.

**Known limitations.** Two folds, one of which contains a volatility event that dominates the
tail behaviour being compared. The IC is an average of daily rank correlations across the
positions live on each date, with no adjustment for the overlap between straddles held at the
same time, so it is a diagnostic rather than a test. The grid varies capacity and loss at a fixed
learning rate. And the fills are at the next close, because the option chain carries end-of-
session quotes only - a constraint of the data described in [`02_labels`](02_labels.ipynb).

**Next**: [`08_tabular_dl`](08_tabular_dl.ipynb) fits a neural network to the same rows. It is
trained by gradient descent on a loss the notebook chooses, so the comparison above can be run
again in a different model class - and if the tails are what separates the objectives here, the
same ordering should appear there.
![notebook output](figures/p1_1.png)
![notebook output](figures/p1_2.png)