Regularized Logistic Regression for ETF Return Direction
Notebook Machine Learning for Trading
Summary
This notebook recasts prediction of 21-session ETF forward returns as a binary task: positive returns are labeled up, and all others down. It fits L2- and L1-regularized logistic regression using chronological walk-forward folds with a purge gap. Scaling is learned within each training fold and then applied to its validation data, reducing leakage. Evaluation includes accuracy, precision, recall, F1, ROC-AUC, log loss, calibration, and coefficient-based feature analysis.
The notebook discusses class imbalance, solver choices, probability calibration, and the potential use of probabilities for position sizing. Its takeaways caution that AUC remains near a constant-score reference, so the exercise does not demonstrate a deployable directional edge. Calibration adjustments may not improve probabilities, and walk-forward folds are for development rather than final performance claims. Any strategy conclusion requires a separate holdout untouched by feature, parameter, and model-selection decisions.
Key ideas
- Convert continuous forward returns into a binary direction target for classification.
- Fit scaling and regularized logistic models separately within each chronological training fold.
- Use ranking and probability metrics alongside accuracy to assess classification quality.
- Check probability calibration before using model probabilities to size positions.
- Reserve an untouched final holdout for performance claims after all development choices are fixed.
Tags
Full text
# Logistic Regression for Return Direction Prediction
# Logistic Regression for Return Direction Prediction
**Docker image**: `ml4t`
**Purpose**: apply logistic regression to predict the direction of 21-day forward
returns on the ETF panel, using the same walk-forward folds as
`02_regularization_paths`.
**Learning objectives**
- Convert a continuous return label into a binary classification target
- Fit L2- and L1-regularized logistic regression on walk-forward folds
- Evaluate with classification metrics: accuracy, AUC-ROC, precision, recall
- Assess probability calibration and its importance for position sizing
- Compare feature importance from L1 logistic coefficients
**Book reference**: Section 11.3 - Predicting Direction with Logistic Regression.
**Prerequisites**
- Ch7 21-day forward return labels at `case_studies/etfs/labels/fwd_ret_21d.parquet`
- Ch8 ETF features at `case_studies/etfs/features/financial.parquet`
- Walk-forward CV configuration in `case_studies/etfs/config/setup.yaml`
**Downstream**: `04_nested_cv_hpo` (HPO with nested CV), `06_conformal_prediction`
(uncertainty quantification), Ch12 (gradient boosting on the same task).
## Setup
```python
"""Logistic Regression for Return Direction Prediction - classify direction and calibration."""
import hashlib
import inspect
from importlib.metadata import version
import joblib
import matplotlib.pyplot as plt
import numpy as np
import polars as pl
import sklearn
from IPython.display import Markdown, display
from matplotlib.patches import Patch
from sklearn.calibration import CalibratedClassifierCV, calibration_curve
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (
accuracy_score,
classification_report,
confusion_matrix,
f1_score,
log_loss,
precision_recall_curve,
precision_score,
recall_score,
roc_auc_score,
roc_curve,
)
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from utils.cv_splits import generate_cv_splits
from utils.modeling import array_sha256, canonical_sha256, file_sha256
from utils.paths import display_path, get_case_study_dir, get_chapter_dir, get_output_dir
from utils.reproducibility import set_global_seeds
from utils.style import COLORS, show_with_alt
```
```python
SEED = 42
MAX_SYMBOLS = 0
RETRAIN = False
```
```python
RANDOM_SEED = SEED
LABEL_HORIZON_SESSIONS = 21
OUTER_LABEL_BUFFER = "21D"
L2_C_VALUES = [0.001, 0.01, 0.1, 1.0, 10.0, 100.0]
L1_C_VALUES = [0.001, 0.01, 0.1, 1.0, 10.0]
MODEL_MAX_ITER = 1000
CACHE_SCHEMA_VERSION = 2
set_global_seeds(SEED)
```
## Load Features and Labels
We reuse the same pre-computed features (Ch8) and 21-day forward return labels
(Ch7) as `02_regularization_paths`. The only difference is that we convert the
continuous return into a binary direction target.
```python
CASE_DIR = get_case_study_dir("etfs")
FEATURES_PATH = CASE_DIR / "features" / "financial.parquet"
LABELS_PATH = CASE_DIR / "labels" / "fwd_ret_21d.parquet"
assert FEATURES_PATH.exists(), (
f"Features not found: {FEATURES_PATH}\nRun the Ch8 ETF features notebook first."
)
assert LABELS_PATH.exists(), (
f"Labels not found: {LABELS_PATH}\nRun the Ch7 ETF labels notebook first."
)
features_df = pl.read_parquet(FEATURES_PATH).with_columns(pl.col("timestamp").cast(pl.Date))
labels_df = pl.read_parquet(LABELS_PATH).with_columns(pl.col("timestamp").cast(pl.Date))
```
### Create Binary Direction Target
We convert the continuous 21-day forward return into a binary label:
1 if the return is positive (up), 0 otherwise (down). This transforms
the regression problem from NB02 into a classification problem.
```python
RETURN_COL = "fwd_ret_21d"
TARGET_COL = "direction"
ASSET_COL = "symbol"
df = features_df.join(labels_df, on=["timestamp", ASSET_COL], how="inner").with_columns(
(pl.col(RETURN_COL) > 0).cast(pl.Int32).alias(TARGET_COL),
)
META_COLS = {"timestamp", ASSET_COL, RETURN_COL, TARGET_COL}
FEATURE_COLS = sorted(c for c in df.columns if c not in META_COLS)
# Drop features that are entirely null (can happen with reduced test universes)
all_null = [c for c in FEATURE_COLS if df[c].null_count() == df.height]
if all_null:
print(f"Dropping {len(all_null)} all-null features: {all_null}")
df = df.drop(all_null)
FEATURE_COLS = [c for c in FEATURE_COLS if c not in all_null]
# Replace inf/NaN with null, then drop all nulls
df = df.with_columns(
[
pl.when(pl.col(c).is_nan() | pl.col(c).is_infinite())
.then(None)
.otherwise(pl.col(c))
.alias(c)
for c in FEATURE_COLS
]
)
df = df.drop_nulls(subset=FEATURE_COLS + [TARGET_COL]).sort(["timestamp", ASSET_COL])
if MAX_SYMBOLS > 0:
assets = df[ASSET_COL].unique().sort().head(MAX_SYMBOLS).to_list()
df = df.filter(pl.col(ASSET_COL).is_in(assets))
print(f"Shape: {df.height:,} rows x {len(FEATURE_COLS)} features")
print(f"Assets: {df[ASSET_COL].n_unique()}")
print(f"Date range: {df['timestamp'].min()} to {df['timestamp'].max()}")
```
### Class Balance
A significant imbalance would require stratified sampling or reweighting.
With 21-day returns on diversified ETFs over a long-term bull market, we
expect a moderate positive skew.
```python
up_count = df.filter(pl.col(TARGET_COL) == 1).height
down_count = df.filter(pl.col(TARGET_COL) == 0).height
print(f"Up (1): {up_count:>8,} ({up_count / df.height:.1%})")
print(f"Down (0): {down_count:>8,} ({down_count / df.height:.1%})")
```
## Walk-Forward Cross-Validation Setup
We reuse the same walk-forward splits from `setup.yaml` as
`02_regularization_paths`: rolling train/validation windows with purge gap.
```python
splits = sorted(
generate_cv_splits(
df,
case_study_id="etfs",
label_buffer=OUTER_LABEL_BUFFER,
date_col="timestamp",
),
key=lambda split: split["val_start"],
)
features_array = df.select(FEATURE_COLS).to_numpy()
target_array = df[TARGET_COL].to_numpy()
return_array = df[RETURN_COL].to_numpy()
dates_np = df["timestamp"].to_numpy()
cv_splits = []
for s in splits:
tr_start, tr_end = np.datetime64(s["train_start"]), np.datetime64(s["train_end"])
te_start, te_end = np.datetime64(s["val_start"]), np.datetime64(s["val_end"])
train_idx = np.where((dates_np >= tr_start) & (dates_np <= tr_end))[0]
test_idx = np.where((dates_np >= te_start) & (dates_np <= te_end))[0]
cv_splits.append((train_idx, test_idx))
train_sizes = [len(tr) for tr, _ in cv_splits]
test_sizes = [len(te) for _, te in cv_splits]
if cv_splits:
print(
f"{len(cv_splits)} walk-forward folds - train size "
f"{min(train_sizes):,}-{max(train_sizes):,}, validation size "
f"{min(test_sizes):,}-{max(test_sizes):,}"
)
else:
print("0 walk-forward folds - every candidate split failed the train/validation size gate")
```
## Helper Functions
### Solver Selection
sklearn's `LogisticRegression` offers several solvers optimized for different
penalty structures:
| Solver | L1 | L2 | ElasticNet | Multinomial | Best for |
|------------|----|----|------------|-------------|----------------------------------------------|
| `lbfgs` | no | yes | no | yes | Default L2; fast quasi-Newton |
| `liblinear` | yes | yes | no | no | Pure L1 binary; coordinate descent |
| `saga` | yes | yes | yes | yes | ElasticNet, large $n$, `sample_weight` |
We use `lbfgs` for L2 (fast, numerically stable), `liblinear` for pure L1
binary classification, and `saga` only when ElasticNet mixing is needed.
```python
def evaluate_classification(y_true: np.ndarray, y_pred: np.ndarray, y_prob: np.ndarray) -> dict:
"""Classification metrics: accuracy, precision, recall, F1, AUC-ROC, log-loss."""
return {
"accuracy": accuracy_score(y_true, y_pred),
"precision": precision_score(y_true, y_pred, zero_division=0),
"recall": recall_score(y_true, y_pred, zero_division=0),
"f1": f1_score(y_true, y_pred, zero_division=0),
# A fold where every label is the same class has no ROC curve and no
# log-loss reference. NaN says that; 0.5 would claim the model was
# measured and found no better than a coin.
"auc_roc": roc_auc_score(y_true, y_prob) if len(np.unique(y_true)) > 1 else np.nan,
"log_loss": log_loss(y_true, y_prob) if len(np.unique(y_true)) > 1 else np.nan,
}
```
### Logistic Estimator
Penalty-specific solver choices live in one constructor so the cache signature
and every fold share the same training contract.
```python
def build_logistic_model(l1_ratio: float, C: float) -> LogisticRegression:
"""Construct the penalty-specific logistic estimator."""
if l1_ratio == 0.0:
penalty, solver = "l2", "lbfgs"
model_kwargs = {}
elif l1_ratio == 1.0:
penalty, solver = "l1", "liblinear"
model_kwargs = {}
else:
penalty, solver = "elasticnet", "saga"
model_kwargs = {"l1_ratio": l1_ratio}
return LogisticRegression(
penalty=penalty,
C=C,
solver=solver,
max_iter=MODEL_MAX_ITER,
random_state=RANDOM_SEED,
**model_kwargs,
)
```
### One Walk-Forward Fold
Every fold learns preprocessing and model state from its training rows only.
```python
def fit_logistic_fold(
train_idx: np.ndarray,
test_idx: np.ndarray,
l1_ratio: float,
C: float,
) -> tuple[dict, np.ndarray, dict, dict]:
"""Fit and evaluate one chronological fold."""
X_tr, X_te = features_array[train_idx], features_array[test_idx]
y_tr, y_te = target_array[train_idx], target_array[test_idx]
scaler = StandardScaler()
X_tr_s = scaler.fit_transform(X_tr)
X_te_s = scaler.transform(X_te)
model = build_logistic_model(l1_ratio, C)
model.fit(X_tr_s, y_tr)
y_pred = model.predict(X_te_s)
y_prob = model.predict_proba(X_te_s)[:, 1]
metrics = evaluate_classification(y_te, y_pred, y_prob)
prediction = {
"y_true": y_te,
"y_pred": y_pred,
"y_prob": y_prob,
"returns": return_array[test_idx],
}
return metrics, model.coef_.ravel().copy(), prediction, {"model": model, "scaler": scaler}
```
### Walk-Forward Logistic CV
The wrapper collects fold metrics, predictions, coefficients, and fitted state.
```python
def cross_validate_logistic(
l1_ratio: float = 0.0,
C: float = 1.0,
) -> tuple[list[dict], list[np.ndarray], list[dict], list[dict]]:
"""Run logistic regression on every walk-forward fold."""
results, coefficients, predictions, fold_models = [], [], [], []
for i, (train_idx, test_idx) in enumerate(cv_splits):
metrics, coefficient, prediction, fold_model = fit_logistic_fold(
train_idx,
test_idx,
l1_ratio,
C,
)
metrics["fold"] = i + 1
results.append(metrics)
coefficients.append(coefficient)
predictions.append(prediction)
fold_models.append(fold_model)
return results, coefficients, predictions, fold_models
```
## Model Cache
Training results are cached to disk so subsequent runs skip model fitting.
Its signature binds the cleaned arrays, split state, model configuration, source
implementation, and relevant library versions. Any semantic change therefore
triggers a genuine retrain rather than silently reusing stale predictions.
Set `RETRAIN = True` to force retraining even when the hashes match.
The three digests the contract is built from are shared with the other notebooks
that cache a fit, in `utils.modeling`: a content hash of each input file, an
array hash that also covers filtering, symbol limits, row order, feature order
and cleaning semantics, and a canonical hash of the nested contract itself.
```python
def assess_cache(cached: dict, expected_signature: dict) -> tuple[list[str], bool]:
"""Return missing result keys and whether the training signature matches."""
required = {"l2_all", "l2_summary", "best_l2_C", "l1_all", "l1_summary", "best_l1_C"}
return sorted(required - set(cached)), cached.get("input_signature") == expected_signature
```
The training contract records every choice that can change fitted predictions.
A canonical digest keeps the cache comparison compact and deterministic.
```python
RESULTS_DIR = get_output_dir(11, "03_logistic_classification")
RESULTS_PATH = RESULTS_DIR / "cv_results.joblib"
SETUP_PATH = CASE_DIR / "config" / "setup.yaml"
cv_indices = np.concatenate(
[np.concatenate([train_idx, [-1], test_idx, [-2]]) for train_idx, test_idx in cv_splits]
).astype(np.int64)
```
The data and split contracts include both declared configuration and the resolved
rows. This catches changes in cleaning, symbol limits, calendar mapping, or fold order.
```python
DATA_CONTRACT = {
"case_study_id": "etfs",
"date_column": "timestamp",
"symbol_column": ASSET_COL,
"return_column": RETURN_COL,
"target_column": TARGET_COL,
"target_rule": "direction = int(fwd_ret_21d > 0)",
"feature_columns": FEATURE_COLS,
"cleaning": "inner join; drop all-null features; finite-only; drop nulls; sort date-symbol",
"scaling": "fold-local StandardScaler fit on training rows only",
"max_symbols": MAX_SYMBOLS,
"symbol_subset": "sorted unique symbols, first max_symbols; zero means all",
}
SPLIT_CONTRACT = {
"source_config_sha256": file_sha256(SETUP_PATH),
"label_buffer": OUTER_LABEL_BUFFER,
"label_horizon_sessions": LABEL_HORIZON_SESSIONS,
"direction": "chronological ascending after sort by val_start",
"resolved_windows": splits,
"resolved_indices_sha256": array_sha256(cv_indices),
}
```
The model contract binds optimization, deterministic selection, implementation
source, and relevant dependency versions.
```python
MODEL_CONTRACT = {
"l2": {"C": L2_C_VALUES, "l1_ratio": 0.0, "penalty": "l2", "solver": "lbfgs"},
"l1": {"C": L1_C_VALUES, "l1_ratio": 1.0, "penalty": "l1", "solver": "liblinear"},
"max_iter": MODEL_MAX_ITER,
"random_seed": RANDOM_SEED,
"selection": {
"metric": "mean fold AUC-ROC",
"aggregation": "unweighted arithmetic mean across folds",
"direction": "maximize",
"tie_break": "smallest C",
},
}
```
```python
TRAINING_CONTRACT = {
"schema_version": CACHE_SCHEMA_VERSION,
"data": DATA_CONTRACT,
"splits": SPLIT_CONTRACT,
"models": MODEL_CONTRACT,
"implementation": {
"notebook_source_sha256": file_sha256(
get_chapter_dir(11) / "03_logistic_classification.py"
),
"splitter_source_sha256": hashlib.sha256(
inspect.getsource(generate_cv_splits).encode()
).hexdigest(),
},
"versions": {
"numpy": np.__version__,
"polars": pl.__version__,
"scikit_learn": sklearn.__version__,
"joblib": joblib.__version__,
"ml4t_diagnostic": version("ml4t-diagnostic"),
"exchange_calendars": version("exchange-calendars"),
"pandas_market_calendars": version("pandas-market-calendars"),
},
}
```
Input hashes bind the contract to the exact cleaned arrays. The complete contract
remains inside the cache so a verifier can inspect what produced each result.
```python
TRAINING_CONTRACT_SHA256 = canonical_sha256(TRAINING_CONTRACT)
INPUT_SIGNATURE = {
"training_contract_sha256": TRAINING_CONTRACT_SHA256,
"features_sha256": file_sha256(FEATURES_PATH),
"labels_sha256": file_sha256(LABELS_PATH),
"features_array_sha256": array_sha256(features_array),
"target_array_sha256": array_sha256(target_array),
"return_array_sha256": array_sha256(return_array),
"dates_sha256": array_sha256(dates_np),
"training_contract": TRAINING_CONTRACT,
}
```
```python
NEED_TRAINING = RETRAIN or not RESULTS_PATH.exists()
if not NEED_TRAINING:
_cached = joblib.load(RESULTS_PATH)
missing_keys, signature_matches = assess_cache(_cached, INPUT_SIGNATURE)
NEED_TRAINING = bool(missing_keys or not signature_matches)
if NEED_TRAINING:
reason = f"missing keys {missing_keys}" if missing_keys else "training signature changed"
print(f"Ignoring stale model cache: {reason}.")
else:
l2_all = _cached["l2_all"]
l2_summary = _cached["l2_summary"]
best_l2_C = _cached["best_l2_C"]
l1_all = _cached["l1_all"]
l1_summary = _cached["l1_summary"]
best_l1_C = _cached["best_l1_C"]
print(f" L2: {len(l2_all)} C values | L1: {len(l1_all)} C values")
del _cached
if NEED_TRAINING:
print("Training models (RETRAIN=True or no cache found)...")
```
## L2 Regularized Logistic Regression (Ridge)
The regularization parameter $C$ is the *inverse* of penalty strength: small
$C$ means a strong penalty and a heavily shrunk model, large $C$ means an
almost unpenalized one. The sweep below spans several orders of magnitude in
both directions, because where the useful range sits is a property of the data
rather than something to carry over from another dataset.
```python
if NEED_TRAINING:
l2_all = {}
for C in L2_C_VALUES:
res, coeffs, preds, models = cross_validate_logistic(l1_ratio=0.0, C=C)
l2_all[C] = {
"results": pl.DataFrame(res),
"coeffs": np.array(coeffs),
"predictions": preds,
"models": models,
}
l2_summary = pl.DataFrame(
[
{
"C": C,
"mean_acc": d["results"]["accuracy"].mean(),
"mean_auc": d["results"]["auc_roc"].mean(),
"mean_f1": d["results"]["f1"].mean(),
}
for C, d in l2_all.items()
]
).sort(["mean_auc", "C"], descending=[True, False])
best_l2_C = l2_summary.row(0, named=True)["C"]
```
```python
print(
f"Best L2 C: {best_l2_C} AUC: {l2_summary.filter(pl.col('C') == best_l2_C)['mean_auc'].item():.4f}"
)
l2_summary
```
### Effect of Class Weighting
`class_weight='balanced'` scales each class's loss contribution inversely
by its frequency. We evaluate the unweighted and balanced objectives on the
latest walk-forward fold. This comparison isolates the effect of weighting
while keeping the training window, regularization, and decision threshold fixed.
```python
train_idx_cw, test_idx_cw = cv_splits[-1]
scaler_cw = StandardScaler()
X_tr_cw = scaler_cw.fit_transform(features_array[train_idx_cw])
X_te_cw = scaler_cw.transform(features_array[test_idx_cw])
y_tr_cw, y_te_cw = target_array[train_idx_cw], target_array[test_idx_cw]
rows_cw = []
for cw_label, cw_val in [("None", None), ("Balanced", "balanced")]:
m = LogisticRegression(
C=best_l2_C,
solver="lbfgs",
max_iter=MODEL_MAX_ITER,
random_state=RANDOM_SEED,
class_weight=cw_val,
)
m.fit(X_tr_cw, y_tr_cw)
yp = m.predict_proba(X_te_cw)[:, 1]
rows_cw.append(
{
"class_weight": cw_label,
**evaluate_classification(y_te_cw, (yp >= 0.5).astype(int), yp),
}
)
cw_df = pl.DataFrame(rows_cw)
cw_df
```
Read the three metrics against each other rather than one at a time. Recall and
F1 describe a decision rule at a fixed threshold, so reweighting the classes
moves them. AUC describes the ranking underneath, which reweighting leaves
largely alone: the same observations are still ordered the same way.
```python
cw_default, cw_balanced = cw_df.iter_rows(named=True)
display(
Markdown(
f"The balanced objective changes majority-class recall from "
f"**{cw_default['recall']:.3f}** to **{cw_balanced['recall']:.3f}** and F1 from "
f"**{cw_default['f1']:.3f}** to **{cw_balanced['f1']:.3f}**, while AUC moves only "
f"from **{cw_default['auc_roc']:.3f}** to **{cw_balanced['auc_roc']:.3f}**. "
"Class weighting therefore changes the fitted decision rule materially even when "
"rank discrimination changes little."
)
)
```
## L1 Regularized Logistic Regression (LASSO)
L1 regularization drives some coefficients to exactly zero, performing
automatic feature selection - the same sparsity effect we saw with LASSO
regression in NB02.
```python
if NEED_TRAINING:
l1_all = {}
for C in L1_C_VALUES:
res, coeffs, preds, models = cross_validate_logistic(l1_ratio=1.0, C=C)
l1_all[C] = {
"results": pl.DataFrame(res),
"coeffs": np.array(coeffs),
"predictions": preds,
"models": models,
}
l1_summary = pl.DataFrame(
[
{
"C": C,
"mean_acc": d["results"]["accuracy"].mean(),
"mean_auc": d["results"]["auc_roc"].mean(),
"mean_f1": d["results"]["f1"].mean(),
"n_nonzero": int((np.abs(d["coeffs"]).mean(axis=0) > 1e-8).sum()),
}
for C, d in l1_all.items()
]
).sort(["mean_auc", "C"], descending=[True, False])
best_l1_C = l1_summary.row(0, named=True)["C"]
```
```python
print(
f"Best L1 C: {best_l1_C} AUC: {l1_summary.filter(pl.col('C') == best_l1_C)['mean_auc'].item():.4f}"
)
l1_summary
```
### Save / Load Cache
```python
if NEED_TRAINING:
RESULTS_DIR.mkdir(parents=True, exist_ok=True)
joblib.dump(
{
"input_signature": INPUT_SIGNATURE,
"l2_all": l2_all,
"l2_summary": l2_summary,
"best_l2_C": best_l2_C,
"l1_all": l1_all,
"l1_summary": l1_summary,
"best_l1_C": best_l1_C,
},
RESULTS_PATH,
)
print(f"Saved results to {display_path(RESULTS_PATH)}")
```
## Model Comparison
We compare the selected L2 and L1 configurations against a naive baseline that
always predicts the majority class of its own training fold.
```python
rows = []
for label, res_df in [
(f"Logistic L2 (C={best_l2_C})", l2_all[best_l2_C]["results"]),
(f"Logistic L1 (C={best_l1_C})", l1_all[best_l1_C]["results"]),
]:
rows.append(
{
"Model": label,
"Accuracy": round(res_df["accuracy"].mean(), 4),
"AUC-ROC": round(res_df["auc_roc"].mean(), 4),
"Precision": round(res_df["precision"].mean(), 4),
"Recall": round(res_df["recall"].mean(), 4),
"F1": round(res_df["f1"].mean(), 4),
"Log-Loss": round(res_df["log_loss"].mean(), 4),
}
)
```
The baseline has to be fitted like a model, fold by fold. Its predicted class
and its predicted probability both come from the training rows of that fold
only, so it never sees the validation window it is scored on.
```python
baseline_folds = []
for fold, (train_idx, test_idx) in enumerate(cv_splits, start=1):
train_rate = float(target_array[train_idx].mean())
majority_class = int(train_rate >= 0.5)
baseline_pred = np.full(len(test_idx), majority_class, dtype=np.int32)
baseline_prob = np.full(len(test_idx), train_rate)
baseline_folds.append(
{
"fold": fold,
**evaluate_classification(target_array[test_idx], baseline_pred, baseline_prob),
}
)
baseline_summary = pl.DataFrame(baseline_folds).select(pl.exclude("fold")).mean().row(0, named=True)
rows.append(
{
"Model": "Naive (majority class)",
"Accuracy": round(baseline_summary["accuracy"], 4),
"AUC-ROC": round(baseline_summary["auc_roc"], 4),
"Precision": round(baseline_summary["precision"], 4),
"Recall": round(baseline_summary["recall"], 4),
"F1": round(baseline_summary["f1"], 4),
"Log-Loss": round(baseline_summary["log_loss"], 4),
}
)
comparison = pl.DataFrame(rows)
comparison
```
The dashed line in each panel is what the majority-class baseline scores on
that panel's metric. Only AUC has a fixed no-information value; for accuracy
and F1 the reference depends on how often the label is Up.
```python
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
model_names = comparison["Model"].to_list()[:2]
x = np.arange(len(model_names))
bar_colors = [COLORS["blue"], COLORS["amber"]]
for ax, metric, title in zip(
axes,
["Accuracy", "AUC-ROC", "F1"],
["Accuracy", "AUC-ROC", "F1 Score"],
strict=False,
):
vals = comparison[metric].to_list()[:2]
bars = ax.barh(x, vals, color=bar_colors)
for bar, val in zip(bars, vals, strict=False):
ax.text(
val + max(vals) * 0.02,
bar.get_y() + bar.get_height() / 2,
f"{val:.3f}",
va="center",
fontsize=9,
)
ax.set_yticks(x)
ax.set_yticklabels(model_names)
ax.set_xlabel(metric)
ax.set_title(title)
baseline_value = comparison.filter(pl.col("Model") == "Naive (majority class)")[metric].item()
ax.axvline(
baseline_value,
ls="--",
color=COLORS["neutral"],
alpha=0.7,
label="Naive baseline",
)
ax.set_xlim(0, max(max(vals), baseline_value) * 1.18)
ax.legend(loc="lower right", fontsize=8)
fig.suptitle("Validation metrics by penalty, against the majority-class baseline")
show_with_alt(
fig,
"Three horizontal bar panels - accuracy, AUC-ROC and F1 - each showing the L2 "
"and L1 models against a dashed line for the majority-class baseline.",
)
```
Two references matter here and they are different numbers. For AUC the
no-information reference is one half, the score of any constant ranking. For
accuracy it is the majority-class rate, which on a panel of ETFs over a long
expansion is well above one half. A model can beat the first and lose to the
second at the same time, which is why accuracy alone is the wrong headline for
a directional model.
```python
l2_row, l1_row, naive_row = comparison.iter_rows(named=True)
display(
Markdown(
f"L2 and L1 reach validation AUCs of **{l2_row['AUC-ROC']:.3f}** and "
f"**{l1_row['AUC-ROC']:.3f}**, respectively, versus **0.500** for a constant "
f"score. Their accuracies of **{l2_row['Accuracy']:.3f}** and "
f"**{l1_row['Accuracy']:.3f}** remain below the majority-class baseline of "
f"**{naive_row['Accuracy']:.3f}**. The features add little directional ranking "
"power, and changing the penalty does not create a meaningfully different operating point."
)
)
```
## Detailed Classification Analysis
We pool the validation predictions of the selected L2 model across all
walk-forward folds to examine the confusion matrix, the ROC curve and the
precision-recall tradeoff. Pooling is legitimate here because every prediction
in it was made out of sample, on the fold that produced it.
```python
best_preds = l2_all[best_l2_C]["predictions"]
all_y_true = np.concatenate([p["y_true"] for p in best_preds])
all_y_prob = np.concatenate([p["y_prob"] for p in best_preds])
all_y_pred = (all_y_prob >= 0.5).astype(int)
```
### Confusion Matrix
```python
cm = confusion_matrix(all_y_true, all_y_pred)
# sklearn returns rows = actual, cols = predicted: cm[i, j] = count(y_true=i, y_pred=j).
pl.DataFrame(
{
"Actual": ["Down", "Up"],
"Predicted Down": [cm[0, 0], cm[1, 0]],
"Predicted Up": [cm[0, 1], cm[1, 1]],
}
)
```
```python
fig, ax = plt.subplots(figsize=(5, 4))
im = ax.imshow(cm, cmap="Blues")
ax.set_xticks([0, 1])
ax.set_yticks([0, 1])
ax.set_xticklabels(["Down", "Up"])
ax.set_yticklabels(["Down", "Up"])
ax.set_xlabel("Predicted")
ax.set_ylabel("Actual")
ax.set_title("Validation confusion matrix at the default threshold")
for i in range(2):
for j in range(2):
ax.text(
j,
i,
f"{cm[i, j]:,}",
ha="center",
va="center",
color="white" if cm[i, j] > cm.max() / 2 else "black",
)
show_with_alt(
fig,
"Two-by-two confusion matrix of actual against predicted direction, each cell "
"labelled with its count and shaded by size.",
)
```
The default threshold of one half is a convention, not a decision. Applied to a model whose
probabilities cluster near the base rate, it sends almost everything to the majority class,
which is why the matrix is lopsided. A threshold is a choice about the cost of each kind of
error, and leaving it at the default is making that choice by not making it.
The confusion matrix counts one particular decision rule: predict Up whenever
the estimated probability clears one half. That threshold is a choice, and it
interacts with how often the label is Up in the first place.
```python
actual_up_rate = all_y_true.mean()
predicted_up_rate = all_y_pred.mean()
display(
Markdown(
f"The model predicts Up on **{predicted_up_rate:.1%}** of validation observations "
f"versus a realized Up rate of **{actual_up_rate:.1%}**. This imbalance explains why "
"false positives outnumber false negatives. Section 11.3 discusses choosing a trading "
"threshold separately from estimating probabilities."
)
)
```
### ROC Curve
The ROC curve plots the true positive rate against the false positive rate as
the decision threshold sweeps from one extreme to the other, so it describes
the ranking rather than any one threshold. The area under it is the probability
that a randomly chosen Up observation is scored above a randomly chosen Down
one; a coin flip scores one half.
```python
auc_val = roc_auc_score(all_y_true, all_y_prob)
fpr, tpr, _ = roc_curve(all_y_true, all_y_prob)
fig, ax = plt.subplots(figsize=(6, 5))
ax.plot(fpr, tpr, linewidth=2, label=f"Logistic L2 (AUC = {auc_val:.3f})")
ax.plot([0, 1], [0, 1], ls="--", color=COLORS["neutral"], label="Random")
ax.set_xlabel("False Positive Rate")
ax.set_ylabel("True Positive Rate")
ax.set_title("ROC curve on the validation folds")
ax.legend(loc="lower right")
show_with_alt(
fig,
"ROC curve of true positive rate against false positive rate, drawn against "
"the diagonal that a random ranking would trace.",
)
```
The curve sits just above the diagonal along its whole length, and the area under it is
barely above a half. On a two-class problem that is a model which orders the validation
cases hardly better than a coin, and every metric downstream of it has to be read in that
light.
### Precision-Recall Curve
For imbalanced classes or when the cost of false positives is high,
the precision-recall curve is more informative than ROC.
```python
prec, rec, _ = precision_recall_curve(all_y_true, all_y_prob)
baseline = all_y_true.mean()
fig, ax = plt.subplots(figsize=(6, 5))
ax.plot(rec, prec, linewidth=2, label="Logistic L2")
ax.axhline(baseline, ls="--", color=COLORS["neutral"], label=f"Baseline ({baseline:.2f})")
ax.set_xlabel("Recall")
ax.set_ylabel("Precision")
ax.set_title("Precision against recall on the validation folds")
ax.legend()
show_with_alt(
fig,
"Precision against recall, with a dashed horizontal line at the share of "
"observations whose realized direction is up.",
)
```
## L1 Feature Importance
The L1 logistic model zeroes out unimportant features. Comparing the
surviving features with the LASSO regression results from NB02 shows
whether the same signals matter for direction prediction.
```python
l1_coeffs = l1_all[best_l1_C]["coeffs"] # (n_folds, n_features)
mean_abs = np.abs(l1_coeffs).mean(axis=0)
mean_signed = l1_coeffs.mean(axis=0)
importance = pl.DataFrame(
{
"feature": FEATURE_COLS,
"mean_abs_coeff": mean_abs,
"mean_signed_coeff": mean_signed,
}
).sort("mean_abs_coeff", descending=True)
nonzero = importance.filter(pl.col("mean_abs_coeff") > 1e-8)
print(f"Non-zero features: {nonzero.height}/{len(FEATURE_COLS)} at C={best_l1_C}")
nonzero.head(15)
```
```python
top15 = importance.head(15)
fig, ax = plt.subplots(figsize=(8, 5))
colors = [
COLORS["blue"] if v >= 0 else COLORS["copper"] for v in top15["mean_signed_coeff"].to_list()
]
ax.barh(range(len(top15)), top15["mean_abs_coeff"].to_list(), color=colors)
ax.set_yticks(range(len(top15)))
ax.set_yticklabels(top15["feature"].to_list())
ax.invert_yaxis()
ax.set_xlabel("Mean |Coefficient|")
ax.set_title("Mean absolute L1 coefficient by feature, colored by sign")
ax.legend(
handles=[
Patch(color=COLORS["blue"], label="Positive mean coefficient"),
Patch(color=COLORS["copper"], label="Negative mean coefficient"),
],
loc="lower right",
)
show_with_alt(
fig,
"Horizontal bars of mean absolute coefficient for the leading features, "
"coloured by the sign of the mean coefficient.",
)
```
## Probability Calibration
A classifier is calibrated when its stated probabilities match observed
frequencies: among the observations it scores at some probability of going up,
that share should actually go up. Discrimination and calibration are separate
properties. A model can rank correctly and still state probabilities that are
uniformly too confident, which matters as soon as the probability is used to
size a position rather than only to sort.
```python
cal_fracs, cal_means = calibration_curve(all_y_true, all_y_prob, n_bins=10)
fig, ax = plt.subplots(figsize=(6, 5))
ax.plot([0, 1], [0, 1], ls="--", color=COLORS["neutral"], label="Perfect calibration")
ax.plot(cal_means, cal_fracs, "o-", linewidth=2, markersize=6, label="Model")
ax.set_xlabel("Mean Predicted Probability")
ax.set_ylabel("Observed Fraction Positive")
ax.set_title("Observed frequency against predicted probability")
ax.legend()
show_with_alt(
fig,
"Reliability diagram: observed fraction of up outcomes against mean predicted "
"probability, in ten bins, against the diagonal of perfect calibration.",
)
```
A calibrated model would track the diagonal. This one is close to flat: across the whole
range of predicted probabilities the observed frequency stays in a narrow band, so moving
the prediction from low to high barely moves what actually happens.
That is the same finding the ROC curve gave, read a different way. A flat calibration curve
and an area near a half are two views of a model whose probabilities do not separate the
classes - and it is worth seeing both, because a model can be well calibrated and
uninformative, or sharp and badly calibrated, and the two plots catch different failures.
**Interpretation**: A curve above the diagonal means the model is
under-confident (actual positive rate exceeds predicted probability),
while a curve below means over-confidence. The curve can still depart from the
diagonal even though logistic loss directly
optimizes probability estimates. The curve, not the model family, decides
whether confidence-based position sizing is defensible.
### Platt Scaling Correction
`CalibratedClassifierCV(method='sigmoid')` fits a logistic function to the
model's raw outputs - Platt scaling. Its internal folds must also respect
chronology. We therefore use expanding date blocks with a 21-session label purge,
then compare the original and corrected curves on the latest outer fold.
The inner calibration splitter keeps complete timestamp groups together. For a
21-session forward label, the last training label must mature strictly before the
first calibration timestamp.
```python
def expanding_calibration_splits(
dates: np.ndarray,
n_splits: int = 3,
gap_sessions: int = LABEL_HORIZON_SESSIONS,
) -> list[tuple[np.ndarray, np.ndarray]]:
"""Build expanding chronological splits for Platt calibration."""
unique_dates = np.unique(dates)
block_size = max(1, len(unique_dates) // 10)
starts = np.linspace(int(0.55 * len(unique_dates)), int(0.85 * len(unique_dates)), n_splits)
splits_out = []
for start in starts.astype(int):
validation_dates = unique_dates[start : min(start + block_size, len(unique_dates))]
train_end = start - gap_sessions
if train_end <= 0 or not len(validation_dates):
continue
train_dates = unique_dates[:train_end]
label_end = unique_dates[train_end - 1 + gap_sessions]
if label_end >= validation_dates[0]:
raise ValueError("Calibration training labels overlap the validation block")
train_local = np.flatnonzero(np.isin(dates, train_dates))
validation_local = np.flatnonzero(np.isin(dates, validation_dates))
splits_out.append((train_local, validation_local))
return splits_out
```
```python
tr_last, te_last = cv_splits[-1]
y_te_platt = target_array[te_last]
calibration_cv = expanding_calibration_splits(dates_np[tr_last])
base_model = make_pipeline(
StandardScaler(),
LogisticRegression(
C=best_l2_C,
solver="lbfgs",
max_iter=MODEL_MAX_ITER,
random_state=RANDOM_SEED,
),
)
cal_model = CalibratedClassifierCV(base_model, method="sigmoid", cv=calibration_cv)
cal_model.fit(features_array[tr_last], target_array[tr_last])
y_prob_platt = cal_model.predict_proba(features_array[te_last])[:, 1]
```
The uncorrected model for the comparison is fitted on exactly the same rows,
with the same fold-local scaler, and differs only in that no sigmoid is applied
to its output.
```python
scaler_platt = StandardScaler()
X_tr_platt = scaler_platt.fit_transform(features_array[tr_last])
X_te_platt = scaler_platt.transform(features_array[te_last])
orig_model = LogisticRegression(
C=best_l2_C,
solver="lbfgs",
max_iter=MODEL_MAX_ITER,
random_state=RANDOM_SEED,
)
orig_model.fit(X_tr_platt, target_array[tr_last])
y_prob_orig = orig_model.predict_proba(X_te_platt)[:, 1]
frac_orig, mean_orig = calibration_curve(y_te_platt, y_prob_orig, n_bins=10)
frac_platt, mean_platt = calibration_curve(y_te_platt, y_prob_platt, n_bins=10)
```
Both curves are drawn on the same fold so the comparison is like for like. A
calibration correction is fitted, which means it can also fit noise: the
question the diagram answers is whether the corrected curve sits closer to the
diagonal, not whether it was supposed to.
```python
fig, ax = plt.subplots(figsize=(6, 5))
ax.plot([0, 1], [0, 1], ls="--", color=COLORS["neutral"], label="Perfect")
ax.plot(mean_orig, frac_orig, "o-", label="Original", markersize=5)
ax.plot(mean_platt, frac_platt, "s-", label="Platt-scaled", markersize=5)
ax.set_xlabel("Mean Predicted Probability")
ax.set_ylabel("Observed Fraction Positive")
ax.set_title("Observed frequency against predicted probability, after Platt scaling")
ax.legend()
show_with_alt(
fig,
"Two reliability curves on the same fold, one from the raw model and one after "
"sigmoid calibration, against the diagonal of perfect calibration.",
)
```
```python
orig_log_loss = log_loss(y_te_platt, y_prob_orig)
platt_log_loss = log_loss(y_te_platt, y_prob_platt)
calibration_direction = "improves" if platt_log_loss < orig_log_loss else "worsens"
display(
Markdown(
f"On the latest outer fold, sigmoid calibration **{calibration_direction}** "
f"log-loss from **{orig_log_loss:.4f}** to **{platt_log_loss:.4f}**. The curve "
"still decides whether the revised probabilities support confidence-based sizing."
)
)
```
### Hit Rate by Confidence
If the classifier is informative, accuracy should increase with confidence,
measured as the distance of the predicted probability from the midpoint. We bin
the pooled validation predictions into confidence quintiles and check whether
accuracy rises across them.
```python
confidence = np.abs(all_y_prob - 0.5)
bin_edges = np.quantile(confidence, np.linspace(0, 1, 6))
bin_edges[0] = -0.001 # include zero
bin_labels = np.clip(np.digitize(confidence, bin_edges) - 1, 0, 4)
hit_rows = []
for b in range(5):
mask = bin_labels == b
if mask.sum() > 0:
acc = (all_y_pred[mask] == all_y_true[mask]).mean()
hit_rows.append(
{
"quintile": b + 1,
"confidence": f"{bin_edges[b]:.3f}-{bin_edges[b + 1]:.3f}",
"n_samples": int(mask.sum()),
"accuracy": round(acc, 4),
}
)
hit_table = pl.DataFrame(hit_rows)
hit_table
```
What to look for: accuracy rising monotonically from the lowest quintile to the
highest is the property that makes distance from one half usable as a strength
signal. Bumps in the middle mean the ordering is only reliable at the extremes.
```python
hit_rows_named = list(hit_table.iter_rows(named=True))
lower_accuracies = [row["accuracy"] for row in hit_rows_named[:4]]
top_hit = hit_rows_named[-1]
top_share = top_hit["n_samples"] / hit_table["n_samples"].sum()
display(
Markdown(
f"Accuracy spans **{min(lower_accuracies):.1%}-{max(lower_accuracies):.1%}** across "
f"the lower four confidence quintiles and reaches **{top_hit['accuracy']:.1%}** in "
f"the highest-confidence quintile. Restricting decisions to that bucket retains "
f"**{top_share:.1%}** of observations. The threshold therefore trades breadth for a "
"higher observed hit rate, but the irregular lower buckets warn against treating raw "
"probability distance as a perfectly ordered strength signal."
)
)
```
### From Probabilities to Trading Signals
The same model produces different portfolios depending on how probabilities
are converted to positions. We demonstrate three common approaches on the
last fold: threshold-based, probability-weighted, and rank-based.
```python
last_prob = l2_all[best_l2_C]["predictions"][-1]["y_prob"]
n_sig = len(last_prob)
# Three conversion methods
thresh_long = last_prob > 0.55
thresh_short = last_prob < 0.45
prob_signal = last_prob - 0.5
ranks = np.argsort(np.argsort(last_prob)) / n_sig
rank_long = ranks >= 0.8
rank_short = ranks < 0.2
print(f"Last fold: {n_sig:,} predictions\n")
print(
f"Threshold (0.45/0.55): {thresh_long.sum():,} long, {thresh_short.sum():,} short, {(~thresh_long & ~thresh_short).sum():,} flat"
)
print(f"Rank-based (20/20): {rank_long.sum():,} long, {rank_short.sum():,} short")
print(f"Prob-weighted: all {n_sig:,} have non-zero weight (mean: {prob_signal.mean():.4f})")
overlap = (thresh_long & rank_long).sum()
print(
f"\nThreshold-long ∩ rank-long overlap: {overlap}/{thresh_long.sum()} ({overlap / max(thresh_long.sum(), 1):.0%})"
)
```
The three methods produce substantially different portfolios from the same
predictions. The choice depends on strategy constraints (position limits,
turnover targets, liquidity). Chapter 17 develops the full framework.
## Multinomial Extension: Ternary Direction
Binary up/down discards information about return magnitude. A multinomial
logistic model with three classes (bottom/middle/top tercile) retains the
ordinal structure and maps naturally to the softmax formulation in Section 11.3.
```python
tr_last, te_last = cv_splits[-1]
y_ret_train = return_array[tr_last]
y_ret_test = return_array[te_last]
tercile_edges = np.quantile(y_ret_train, [1 / 3, 2 / 3])
y_tern_train = np.digitize(y_ret_train, tercile_edges) # 0=bottom, 1=mid, 2=top
y_tern_test = np.digitize(y_ret_test, tercile_edges)
scaler_mn = StandardScaler()
X_tr_mn = scaler_mn.fit_transform(features_array[tr_last])
X_te_mn = scaler_mn.transform(features_array[te_last])
model_mn = LogisticRegression(
solver="lbfgs",
C=best_l2_C,
max_iter=MODEL_MAX_ITER,
random_state=RANDOM_SEED,
)
model_mn.fit(X_tr_mn, y_tern_train)
y_pred_mn = model_mn.predict(X_te_mn)
ternary_report = classification_report(
y_tern_test,
y_pred_mn,
target_names=["Bottom", "Middle", "Top"],
output_dict=True,
zero_division=0,
)
pl.DataFrame(
[
{"class": name, **ternary_report[name]}
for name in ["Bottom", "Middle", "Top", "macro avg", "weighted avg"]
]
)
```
Compare the per-class recalls rather than the headline accuracy. Equal recall
across the three terciles would mean the model finds extremes as readily as it
finds the middle; unequal recall means it is systematically better at one end
of the return distribution than the other.
```python
top_metrics = ternary_report["Top"]
bottom_metrics = ternary_report["Bottom"]
display(
Markdown(
f"Latest-fold ternary accuracy is **{ternary_report['accuracy']:.1%}** versus a "
f"one-third random reference. Top-tercile recall is **{top_metrics['recall']:.1%}** "
f"and Bottom-tercile recall is **{bottom_metrics['recall']:.1%}**. The asymmetric "
"errors show that a vanilla multinomial loss does not recover return extremes evenly; "
"class reweighting, an ordinal objective, or a cross-sectional rank mapping would encode "
"that trading objective more directly."
)
)
```
## Key Takeaways
1. **Direction prediction** reframes the return forecasting problem as binary
classification. The majority-class baseline remains a demanding accuracy
reference because positive 21-day ETF returns are more common.
2. **L2 logistic regression** keeps all features and is well-suited when many
correlated signals each contribute a small amount. L1 can create sparsity
under stronger regularization, but the AUC-optimal fit retains most features.
3. **AUC-ROC** is more informative than raw accuracy because it evaluates
discrimination across all probability thresholds. Here it remains close to
the constant-score reference, so the exercise demonstrates workflow rather
than a deployable directional edge.
4. **Probability calibration** matters for position sizing. The reliability
diagrams show gaps between predicted probabilities and observed frequencies.
Chronological Platt scaling reshapes those probabilities but does not guarantee
improvement, so evaluate calibration on data that follows every fitted step.
5. **Walk-forward validation is not a final holdout.** These folds support
model comparison and diagnostics. A claim about what a strategy would have
earned needs a stretch of data that was held back from every choice made
along the way, including the choice of $C$.
**Next**: `04_nested_cv_hpo` adds Optuna-based hyperparameter optimization
with proper nested cross-validation to control selection bias.
Shown in full with attribution under the source's licence. Licence: MIT
This summary was written by Stratmill's research agent from the original; it is not a copy of the source.