Regularized regressions: Ridge, Lasso, Elastic Net and logistic.

Why penalize a linear model

A cross-sectional linear model predicts the following week's return of each stock as a weighted sum of its features. It is the simplest model you can write and, on financial data, the hardest to beat: the signal is small, almost linear in the ranks of the features and changes over time, and a flexible model spends its flexibility fitting noise. That is why regularized regressions are the baseline of the whole module: every model in lessons 5.5.2 to 5.7.5 has to beat them out of sample to justify its cost.

With $n$ observations, the feature matrix $X$ (standardized) and the label $y$, the five variants of this lesson minimize, as scikit-learn 1.9 writes them:

where $\rho$ is l1_ratio (0.5 by default in ElasticNet). Logistic regression models the probability that the stock ends up above the median, $p = 1/(1 + e^{-x'w})$, minimizes the log-loss and is penalized with C (the inverse of the strength) and l1_ratio.

Two scikit-learn details worth knowing before touching α:

• Ridge does not divide the loss by $n$; Lasso and Elastic Net do. With orthogonal features of variance 1, Ridge multiplies each OLS coefficient by $1/(1 + \alpha/n)$: with 145,600 rows, alpha=1e5 leaves it at 59%. Lasso applies soft thresholding: it subtracts α from its absolute value and sets it to zero if it does not reach it, $w_j = \operatorname{sign}(\hat w_j)\max(|\hat w_j| - \alpha, 0)$. If the label is in basis points, so is the Lasso's α: alpha=1 takes 1 bp off each coefficient.
• In logistic regression, penalty has been deprecated since 1.8 and disappears in 1.10: the penalty is chosen with l1_ratio (0 = L2, 1 = L1) and C=np.inf removes the penalty. Old code with penalty="l1" emits a FutureWarning.

Geometry explains the difference between the two penalties. Ridge's is a round ball: the optimum touches the ball at any point, with all the coefficients non-zero and smaller. The Lasso's is a diamond with corners on the axes: the optimum tends to land on a corner, where several coefficients are exactly zero. That is why Ridge shrinks and Lasso selects. Elastic Net has corners (it selects) and curved sides (it spreads the weight among correlated features instead of keeping one at random).

flowchart TB
  A[Features] --> B{Highly correlated}
  B -->|yes| C{Need to choose}
  B -->|no| D[OLS or mild Ridge]
  C -->|no| E[Ridge]
  C -->|yes| F[Elastic Net]
  A --> G{Many candidates}
  G -->|yes| H[Lasso or ENet]

None of the three replaces what really decides the result: features known at their date (lesson 5.2.2), leakage-free normalization (5.2.4) and a validation that respects time (5.4.2 and 5.4.3).

The test bench: the panel with published truth

The data are synthetic: 400 fictitious stocks, 730 weeks from 2012 to 2025 and 21 weekly features, generated with a model whose effects are published in public/curso/datos/nivel-5/README.md. The expected excess return of the following week (alfa_verdadero, in panel5_verdad.parquet) depends on six effects with coefficients of the order of 5 bp per week: momentum, one-week reversal (stronger with low volume), value, earnings surprise that decays, low volatility and sentiment. Another 14 features have no effect of their own, five of them pure noise. The leak revision_analistas is left out (exercise 1 uses it). Nothing that comes out here is evidence about real markets; it is a laboratory where you know what each method should find.

The features are transformed into cross-sectional normal ranks each week (lesson 5.2.4), with 0 when the data point is missing (sentimiento is missing in 32% of the rows), and the label is the following week's return minus the mean of its week, in basis points.

import numpy as np
import pandas as pd
from scipy.stats import norm

f = pd.read_parquet("panel5_features.parquet")
e = pd.read_parquet("panel5_etiquetas.parquet")
v = pd.read_parquet("panel5_verdad.parquet")
d = (f.merge(e[["fecha", "ticker", "ret_sig_semana"]],
             on=["fecha", "ticker"])
      .merge(v, on=["fecha", "ticker"])
      .dropna(subset=["ret_sig_semana"])
      .reset_index(drop=True))
FE = [c for c in f.columns[3:] if c != "revision_analistas"]

def g(s):
    """Normal rank within each week (lesson 5.2.4); missing, 0."""
    por = s.groupby(d["fecha"])
    u = (por.rank() - 0.5) / por.transform("count")
    return pd.Series(norm.ppf(u), index=s.index).fillna(0.0)

def sin_media(s):
    return s - s.groupby(d["fecha"]).transform("mean")

G = pd.DataFrame({k: g(d[k]) for k in FE})
y = 1e4 * sin_media(d["ret_sig_semana"])     # bp per week
sem = d["fecha"].to_numpy()
ent = sem < np.datetime64("2019-01-01")       # 2012-2018
print(G.shape, "training rows:", ent.sum())
print(f"standard deviation of y: {y[ent].std():.0f} bp")

(291600, 21) training rows: 145600
standard deviation of y: 365 bp

The label has a standard deviation of 365 bp per week and the effects are about 5 bp per unit of feature: a signal-to-noise ratio of one to seventy. To know what a linear model should recover in this representation, the next cell projects the true alpha onto the 21 features. Those coefficients are the linear truth: the best a linear model can do on the ranks of these features.

from sklearn.linear_model import LinearRegression

alfa = 1e4 * sin_media(d["alfa_verdadero"])
proy = LinearRegression().fit(G[ent], alfa[ent])
verdad = pd.Series(proy.coef_, index=FE)
print(f"R² of the alpha on G: {proy.score(G[ent], alfa[ent]):.2f}")
print(f"R² of the oracle on y: "
      f"{alfa[ent].var() / y[ent].var():.4f}")
verdad[verdad.abs() > 0.3].round(2)

R² of the alpha on G: 0.61
R² of the oracle on y: 0.0010
ret_1s         -5.13
ret_1m          0.44
mom_12_1        4.65
vol_63d        -3.73
dist_max_52s   -0.49
valor_bp        3.59
ey              0.41
roe            -0.36
sue             2.92
sentimiento     2.93
dtype: float64

Two figures to keep in mind throughout the lesson. The first: a linear model on the ranks explains 61% of the variance of the true alpha; the rest is the interaction of reversal with volume, the decay of the SUE, the saturation of sentiment, the customer-graph component (lesson 5.7.4) and the weekly variation of the coefficients. The second: the oracle, which knows the exact alpha, explains 0.10% of the variance of the label. An out-of-sample R² of 1% on this panel would not be a great model: it would be a leak.

OLS against the truth

With 145,600 rows, OLS has a lot of information. The standard error has to be computed properly: the 400 stocks of the same week share the shocks of their sector, so the rows are not independent. The cell clusters the errors by week (cluster standard errors, lesson 2.1.1).

import statsmodels.api as sm

mco = sm.OLS(y[ent].to_numpy(), sm.add_constant(G[ent])).fit(
    cov_type="cluster", cov_kwds={"groups": sem[ent]})
tabla = pd.DataFrame({"truth": verdad, "ols": mco.params[FE],
                      "t": mco.tvalues[FE]})
tabla.round(2)

                    truth   ols     t
ret_1s              -5.13 -5.10 -2.25
ret_1m               0.44  2.36  0.80
ret_3m               0.16 -1.99 -0.82
mom_12_1             4.65  5.29  2.23
vol_21d             -0.01 -1.26 -0.52
vol_63d             -3.73  1.76  0.69
volumen_rel         -0.00 -0.20 -0.20
rsi_14               0.02 -2.12 -0.85
dist_max_52s        -0.49 -0.65 -0.26
log_capitalizacion   0.04 -0.45 -0.43
valor_bp             3.59  5.90  2.77
ey                   0.41 -4.34 -1.41
roe                 -0.36  2.42  1.09
sue                  2.92  2.04  1.74
dias_desde_informe   0.02  0.23  0.23
sentimiento          2.93  1.91  1.64
ruido_1             -0.00 -1.25 -1.26
ruido_2             -0.00  1.27  1.40
ruido_3              0.07 -0.15 -0.17
ruido_4             -0.00 -0.02 -0.02
ruido_5             -0.05 -0.24 -0.26

The figure sums up the table. Three readings:

• The large effects are recovered. ret_1s comes out at −5.10 with a true −5.13 and mom_12_1 at 5.29 with 4.65. But with clustered errors the t of reversal is −2.25, not −5 as the naive formula would give: seven years of data give less independent information than 145,600 rows suggest.
• Collinearity scrambles the coefficients. ey comes out at −4.34 (truth 0.41), valor_bp at 5.90 (truth 3.59) and roe at 2.42 (truth −0.36). The three offset each other: raising one and lowering another barely changes the fit.
• Confounding is not fixed with statistics. vol_63d has a true effect of −3.73 and OLS estimates it at +1.76. Volatile stocks have more beta and the synthetic market rises on average, so the label (which only removes the week's mean, not beta) rewards volatility. It is the same effect the data README documents with Fama-MacBeth. No penalty corrects it; it takes another label (the return net of beta) or another feature.

print(G.loc[ent, ["valor_bp", "ey", "roe", "sue"]].corr().round(2))
print(G.loc[ent, ["ret_1s", "rsi_14", "mom_12_1", "ruido_4"]]
      .corr().round(2))

          valor_bp    ey   roe   sue
valor_bp      1.00  0.71  0.16  0.02
ey            0.71  1.00  0.72  0.18
roe           0.16  0.72  1.00  0.29
sue           0.02  0.18  0.29  1.00
          ret_1s  rsi_14  mom_12_1  ruido_4
ret_1s      1.00    0.55      0.01     0.00
rsi_14      0.55    1.00      0.02     0.01
mom_12_1    0.01    0.02      1.00     0.60
ruido_4     0.00    0.01      0.60     1.00

ey (earnings over market cap) is almost the product of valor_bp (book equity over market cap) and roe (earnings over book equity): correlations of 0.71 and 0.72. rsi_14 is a version of reversal (0.55 with ret_1s) and ruido_4 was built as 0.6 × momentum plus noise (0.60). These are the groups where an unpenalized linear model spreads the weight badly.

Ridge and Lasso: shrink or select

The next cell fits Ridge with $\alpha = 0.5\,n$ (each coefficient multiplied by about two thirds if the features were orthogonal) and Lasso with $\alpha = 1$ (a threshold of 1 bp).

from sklearn.linear_model import Lasso, Ridge

n = ent.sum()
X, Y = G[ent].to_numpy(), y[ent].to_numpy()
tabla = tabla[["truth", "ols"]].copy()
tabla["ridge"] = Ridge(alpha=0.5 * n).fit(X, Y).coef_
tabla["lasso"] = Lasso(alpha=1.0).fit(X, Y).coef_
tabla.round(2)

                    truth   ols  ridge  lasso
ret_1s              -5.13 -5.10  -3.25  -4.41
ret_1m               0.44  2.36  -0.05  -0.00
ret_3m               0.16 -1.99  -0.76  -0.09
mom_12_1             4.65  5.29   2.40   2.92
vol_21d             -0.01 -1.26  -0.13   0.00
vol_63d             -3.73  1.76   0.81   0.12
volumen_rel         -0.00 -0.20  -0.17  -0.00
rsi_14               0.02 -2.12  -1.42  -0.78
dist_max_52s        -0.49 -0.65   0.09  -0.00
log_capitalizacion   0.04 -0.45  -0.28  -0.00
valor_bp             3.59  5.90   1.81   2.12
ey                   0.41 -4.34   0.19   0.00
roe                 -0.36  2.42   0.24   0.00
sue                  2.92  2.04   1.54   1.27
dias_desde_informe   0.02  0.23   0.19   0.00
sentimiento          2.93  1.91   1.10   0.42
ruido_1             -0.00 -1.25  -0.83  -0.24
ruido_2             -0.00  1.27   0.86   0.28
ruido_3              0.07 -0.15  -0.13  -0.00
ruido_4             -0.00 -0.02   0.63   0.00
ruido_5             -0.05 -0.24  -0.15  -0.00

Ridge undoes the offsetting within the value group: ey goes from −4.34 to 0.19 and roe from 2.42 to 0.24, close to their true values. In exchange, everything shrinks (reversal drops to −3.25) and the grouping effect appears: ruido_4, which OLS left at −0.02, receives 0.63 because Ridge spreads the momentum weight between the two correlated features (mom_12_1 falls to 2.40). For a ranking it does not matter that everything shrinks; what matters are the proportions.

Lasso sets 11 of the 21 features to zero, among them ey and roe. But it does not get the list right: it keeps ruido_1 (−0.24) and ruido_2 (0.28), pure noise; rsi_14 (−0.78), which only copies reversal; and vol_63d with the wrong sign (0.12). And it shrinks sentimiento a lot (0.42 against a truth of 2.93).

from sklearn.linear_model import lasso_path

alfas = np.geomspace(6, 0.05, 40)
_, camino, _ = lasso_path(X - X.mean(0), Y - Y.mean(), alphas=alfas)
camino = pd.DataFrame(camino.T, index=alfas.round(3), columns=FE)
entra = (camino.abs() > 0).idxmax()     # first α, from largest to smallest
entra[camino.abs().max() > 0].sort_values(ascending=False).head(12)

ret_1s                5.307
mom_12_1              3.248
sue                   2.541
valor_bp              2.541
rsi_14                1.988
ruido_2               1.216
ruido_1               1.216
sentimiento           1.216
vol_63d               1.076
ret_3m                0.952
ret_1m                0.356
log_capitalizacion    0.356
dtype: float64

The regularization path runs α from largest to smallest and records when each feature enters. The figure at the top draws seven of those paths. The strongest effects enter first (ret_1s at α = 5.3; mom_12_1 at 3.2; sue and valor_bp at 2.5), then rsi_14 (2.0), which has no effect of its own but resembles reversal, and at α = 1.2 ruido_1, ruido_2 and sentimiento enter together: the Lasso does not distinguish a true 3 bp effect from a chance correlation in seven years of data. vol_63d enters afterwards, with the sign of beta.

The Lasso selects noise

That a feature enters the model does not say it has an effect. A cheap way to check is to repeat the selection in different periods and see whether it repeats and with the same sign. The cell fits the same Lasso (α = 1.5) in seven two-year windows, from 2012-2013 to 2024-2025.

anio = pd.DatetimeIndex(sem).year.to_numpy()
signos = {}
for a0 in range(2012, 2026, 2):
    m = (anio >= a0) & (anio <= a0 + 1)
    las = Lasso(alpha=1.5).fit(G[m], y[m])
    signos[a0] = pd.Series(np.sign(las.coef_), index=FE)
sg = pd.DataFrame(signos)
cuenta = pd.DataFrame({"selected": (sg != 0).sum(axis=1),
                       "positive": (sg > 0).sum(axis=1),
                       "negative": (sg < 0).sum(axis=1),
                       "effect": verdad.round(1)})
cuenta.sort_values("selected", ascending=False).head(12)

                    selected  positive  negative  effect
ret_1s                     6         0         6    -5.1
valor_bp                   6         5         1     3.6
log_capitalizacion         5         3         2     0.0
sue                        5         5         0     2.9
ret_3m                     4         1         3     0.2
mom_12_1                   4         4         0     4.7
vol_21d                    4         3         1    -0.0
ruido_1                    4         1         3    -0.0
ret_1m                     3         2         1     0.4
sentimiento                3         3         0     2.9
dias_desde_informe         3         1         2     0.0
vol_63d                    3         3         0    -3.7

Counting how many times it enters is not enough: log_capitalizacion, with no effect, enters in 5 of the 7 windows, more than mom_12_1 (4) or sentimiento (3). What separates them is the sign: the true effects keep it (ret_1s, 6 negative out of 6; sue, 5 positive out of 5; mom_12_1, 4 of 4; sentimiento, 3 of 3) and the features with no effect change it (log_capitalizacion, 3 and 2; ruido_1, 1 and 3; dias_desde_informe, 1 and 2). The exception is the usual one: vol_63d is stable (3 of 3)... with the wrong sign, because the confounding with beta is also stable. Stability rules out noise; it does not rule out a bias.

With 21 features the problem is small. With many candidates, which is the norm when features are generated in bulk, it is serious. The next cell adds the 210 pairwise products of features (231 columns) and looks at what enters first. The only true interaction is ret_1s*volumen_rel: reversal is stronger with low volume.

from itertools import combinations

pares = list(combinations(FE, 2))
nombres = FE + [f"{a}*{b}" for a, b in pares]
X231 = np.column_stack(
    [G.to_numpy()] + [(G[a] * G[b]).to_numpy() for a, b in pares]
).astype(np.float32)
Xc = X231[ent] - X231[ent].mean(0)
al, cam, _ = lasso_path(Xc, Y - Y.mean(), precompute=True,
                        alphas=np.geomspace(6, 1.5, 12))
cam = pd.DataFrame(cam.T, index=al.round(2), columns=nombres)
orden = (cam.abs() > 0).idxmax()[cam.abs().max() > 0]
orden.sort_values(ascending=False).head(14)

ret_1s                        5.29
mom_12_1                      3.62
valor_bp                      2.48
sue                           2.48
rsi_14                        2.19
ret_1m*volumen_rel            2.19
mom_12_1*ey                   2.19
roe*dias_desde_informe        2.19
vol_63d*dias_desde_informe    2.19
vol_63d*ruido_1               2.19
valor_bp*ruido_5              2.19
sue*ruido_1                   1.93
ret_1s*volumen_rel            1.93
ret_3m*mom_12_1               1.93
dtype: float64

At α = 2.19 six false interactions enter at once (among them vol_63d*ruido_1 and valor_bp*ruido_5, noise by construction) and the true one, ret_1s*volumen_rel, enters one step later, together with two other false ones. With 210 candidates, the maximum of 210 chance correlations is as large as a small true effect: it is the multiple testing problem of lesson 2.1.8 in another disguise. float32 in X231 is not a detail: the matrix takes 270 MB instead of 540, and precompute=True makes the Lasso work with $X'X$ (231 × 231) instead of going through the rows at each iteration.

Does it matter out of sample? The cell trains on 2012-2018 (without the last week, whose label ends in 2019: the minimal purge of lesson 5.4.2) and measures the weekly IC from 2019 to 2025.

def ic(score, mask):
    """Rank IC of each week between score and the label."""
    t = pd.DataFrame({"s": score, "f": sem[mask],
                      "y": d.loc[mask, "ret_sig_semana"].to_numpy()})
    t[["s", "y"]] = t.groupby("f")[["s", "y"]].rank()
    return t.groupby("f")[["s", "y"]].corr().xs("s", level=1)["y"]

def resumen(serie):
    se = serie.std() / np.sqrt(len(serie))
    return f"IC {serie.mean():.4f} ± {se:.4f}  t {serie.mean()/se:4.1f}"

purga = sem < np.datetime64("2018-12-24")    # without the last week
pru = sem >= np.datetime64("2019-01-01")
Ytr = y[purga].to_numpy()
for nombre, cols, mod in [
        ("ols 21", slice(0, 21), LinearRegression()),
        ("ols 231", slice(None), LinearRegression()),
        ("ridge 231", slice(None), Ridge(alpha=2.0 * n)),
        ("lasso 231", slice(None), Lasso(alpha=2.0, precompute=True))]:
    mod.fit(X231[purga][:, cols], Ytr)
    s = mod.predict(X231[pru][:, cols])
    print(f"{nombre:10s}", resumen(ic(s, pru)))
print(f"{'oracle':10s}", resumen(ic(alfa[pru].to_numpy(), pru)))

ols 21     IC 0.0195 ± 0.0045  t  4.4
ols 231    IC 0.0133 ± 0.0035  t  3.9
ridge 231  IC 0.0134 ± 0.0036  t  3.7
lasso 231  IC 0.0188 ± 0.0051  t  3.7
oracle     IC 0.0256 ± 0.0052  t  4.9

Adding 210 candidates to OLS costs it a third of the IC (from 0.0195 to 0.0133). Ridge does not fix it (0.0134): it shrinks the 210 noise interactions but does not remove them, and their sum keeps moving the ranking. Lasso does recover almost all of it (0.0188) because it sets most of them to zero. None beats the small model. Regularization is insurance against a bad set of candidates, not a way of turning them into signal; the cheapest insurance is not to generate 210 candidates without a hypothesis behind them (lesson 2.3.5).

Honest validation: walk-forward with α chosen without looking at the test

Up to here α has been set by hand. In practice it is chosen with validation, and the validation has to respect time. The scheme of the next cell, for each year from 2015 to 2025:

flowchart TB
  A[Data up to t−1] --> B[Sub-training]
  A --> C[Validate t−2 and t−1]
  B --> D[IC for each α]
  C --> D
  D --> E[Best α]
  E --> F[Retrain on everything]
  F --> G[Predict year t]

• Training: all the weeks before the year, except the last one (its label ends inside the test year).
• Validation: the last two years of the training set; sub-training ends one week before they begin.
• Selection: the α with the highest mean IC in validation; then the model is retrained on the whole training set and the year is predicted.

It is a walk-forward with an expanding window (lesson 5.4.3) with a purged validation inside (lesson 5.4.2). The label covers one week, so purging one week is enough; with 21-session labels you would have to purge four.

from sklearn.linear_model import ElasticNet

modelos = {
    "ols": (lambda a: LinearRegression(), [None]),
    "ridge": (lambda a: Ridge(alpha=a), [1e3, 1e4, 1e5, 1e6]),
    "lasso": (lambda a: Lasso(alpha=a), [0.1, 0.5, 1, 2, 3]),
    "enet": (lambda a: ElasticNet(alpha=a), [0.2, 1, 2, 4, 6]),
}
semana = np.timedelta64(7, "D")
Xg, yg = G.to_numpy(), y.to_numpy()
pred = {k: np.full(len(d), np.nan) for k in modelos}
alfa_elegido = {k: [] for k in modelos}
for a in range(2015, 2026):
    ini, ini_val = (np.datetime64(f"{a}-01-01"),
                    np.datetime64(f"{a - 2}-01-01"))
    tr = sem < ini - semana                  # one-week purge
    sub = sem < ini_val - semana
    val = tr & (sem >= ini_val)              # the last two years
    te = anio == a
    for k, (hacer, rejilla) in modelos.items():
        mejor = rejilla[0]
        if len(rejilla) > 1:
            notas = [ic(hacer(r).fit(Xg[sub], yg[sub])
                        .predict(Xg[val]), val).mean()
                     for r in rejilla]
            mejor = rejilla[int(np.argmax(notas))]
        alfa_elegido[k].append(mejor)
        pred[k][te] = hacer(mejor).fit(Xg[tr], yg[tr]).predict(Xg[te])
wf = anio >= 2015
ics = {k: ic(p[wf], wf) for k, p in pred.items()}
ics["oracle"] = ic(alfa[wf].to_numpy(), wf)
for k, s in ics.items():
    print(f"{k:8s}", resumen(s))

ols      IC 0.0184 ± 0.0038  t  4.9
ridge    IC 0.0175 ± 0.0039  t  4.5
lasso    IC 0.0185 ± 0.0039  t  4.7
enet     IC 0.0175 ± 0.0040  t  4.4
oracle   IC 0.0269 ± 0.0040  t  6.8

Over 574 test weeks, OLS gives an IC of 0.0184 ± 0.0038, the same as the reference ridge in the data README (0.0185). Ridge (0.0175), Lasso (0.0185) and Elastic Net (0.0175) are within one standard error. The oracle gives 0.0269: the linear model captures about two thirds of the signal that exists, and what is missing is not recovered by penalizing. The table of chosen α values explains why validation does not help.

pd.DataFrame(alfa_elegido, index=range(2015, 2026)).drop(columns="ols")

          ridge  lasso  enet
2015   100000.0    3.0   4.0
2016  1000000.0    3.0   4.0
2017  1000000.0    1.0   2.0
2018    10000.0    0.5   0.2
2019     1000.0    0.1   0.2
2020     1000.0    0.1   0.2
2021    10000.0    0.5   0.2
2022  1000000.0    1.0   2.0
2023  1000000.0    3.0   6.0
2024     1000.0    0.1   0.2
2025     1000.0    0.1   0.2

Ridge's α jumps from the minimum to the maximum of the grid from one year to the next (1,000 in 2019-2020, 1,000,000 in 2022-2023). It is not instability in the method: the validation curve is almost flat. With $n \gg p$, the penalty barely changes the direction of the coefficient vector, and a rank IC depends only on the direction (multiplying all the predictions by 0.5 does not change any ranking). What validation picks is, above all, two years of noise. The fair comparison between models is paired: the same weeks, the week-by-week IC difference and its standard error.

for k in ["ridge", "lasso", "enet", "oracle"]:
    dif = ics[k] - ics["ols"]
    se = dif.std() / np.sqrt(len(dif))
    print(f"{k:8s} {dif.mean():+.4f} ± {se:.4f}  "
          f"t {dif.mean() / se:+.1f}")

ridge    -0.0010 ± 0.0010  t -1.0
lasso    +0.0000 ± 0.0011  t +0.0
enet     -0.0009 ± 0.0012  t -0.8
oracle   +0.0085 ± 0.0042  t +2.0

Pairing reduces the standard error from 0.0038 to 0.0010: the models get the same weeks wrong. Even so, no difference is significant; the oracle's barely is (t = 2.0), which gives an idea of how hard it is to prove out of sample even an improvement that exists. With five models tested, moreover, the best IC is biased upward; the correction for the number of tests (Deflated Sharpe) is lesson 5.4.4.

Logistic regression: direction instead of size

Logistic regression changes the question: instead of the return, the probability that the stock ends up above the median of its week. It has two advantages on financial data: the binary label ignores the tails (a +40% counts the same as a +1%) and the output is a probability that can be calibrated. The cost is throwing away the size information.

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import log_loss

z = (d["ret_sig_semana"]
     > d.groupby("fecha")["ret_sig_semana"].transform("median"))
z = z.astype(int).to_numpy()
p = np.full(len(d), np.nan)
for a in range(2015, 2026):
    tr = sem < np.datetime64(f"{a}-01-01") - semana
    te = anio == a
    lr = LogisticRegression(C=1.0, l1_ratio=0.0, max_iter=500)
    p[te] = lr.fit(Xg[tr], z[tr]).predict_proba(Xg[te])[:, 1]
acierto = ((p[wf] > 0.5) == z[wf]).mean()
print(f"hit rate {100 * acierto:.2f}%  "
      f"log-loss {log_loss(z[wf], p[wf]):.5f} (chance {np.log(2):.5f})")
print("percentiles 1, 50 and 99 of p:",
      np.percentile(p[wf], [1, 50, 99]).round(3))
print(resumen(ic(p[wf], wf)))

hit rate 50.60%  log-loss 0.69314 (chance 0.69315)
percentiles 1, 50 and 99 of p: [0.471 0.5   0.527]
IC 0.0160 ± 0.0036  t  4.4

The hit rate is 50.60%; that of the oracle over the same period, according to the README, is 50.89%. The log-loss, 0.69314, is one hundred-thousandth below that of a coin (ln 2 = 0.69315), and 98% of the probabilities lie between 0.471 and 0.527. That is what a good model looks like when the true R² is 0.1%: almost indistinguishable from chance stock by stock, useful when 400 stocks are ranked at once every week for years. The IC, 0.0160, is somewhat below that of OLS (0.0184): the binary label throws away information.

C=1.0 with 150,000 rows is an almost negligible penalty: in scikit-learn's logistic regression the loss is a sum, not a mean, so the weight of the penalty falls with $n$. With l1_ratio=1.0 (L1) you need solver="saga" or "liblinear"; both are slower at this size.

lr = LogisticRegression(C=1.0, l1_ratio=0.0, max_iter=500)
lr.fit(Xg[ent], z[ent])
comp = pd.DataFrame({"ols_bp": tabla["ols"],
                     "logit_x1000": 1000 * lr.coef_[0]}, index=FE)
print(f"correlation: {comp.corr().iloc[0, 1]:.3f}")
comp.loc[["ret_1s", "mom_12_1", "valor_bp", "sue", "ey"]].round(2)

correlation: 0.933
          ols_bp  logit_x1000
ret_1s     -5.10       -22.33
mom_12_1    5.29        22.88
valor_bp    5.90        30.85
sue         2.04         5.71
ey         -4.34       -14.87

The logistic and OLS coefficients have a correlation of 0.933: the two models see the same thing. The differences are where the tail matters: sue weighs relatively less in the logistic regression (5.71 against 22.33 for reversal, when in OLS it is 2.04 against 5.10), because a good part of its effect arrives in large jumps that the binary label treats like any other rise.

Common mistakes

• Penalizing without standardizing. Ridge and Lasso penalize all coefficients equally; if one feature is in millions and another in fractions of one, the penalty punishes the second. Use normal ranks or StandardScaler inside the pipeline (lesson 5.1.2), fitted only on the training set.
• Reading Ridge's α as the Lasso's. In scikit-learn, one scales with $n$ and the other does not. An alpha=1 in Ridge with 145,600 rows does almost nothing; in Lasso it takes 1 bp off each coefficient.
• Interpreting the Lasso's selection as evidence. Whatever correlates with the residual in that sample enters, noise included. You have to repeat it in other periods and require the same sign.
• Choosing α with shuffled k-fold. In this panel the labels are weekly and do not overlap, but with 21-session labels random k-fold promises 60% hit rates that do not exist (lesson 5.4.1). Always use purged validation.
• Expecting the penalty to correct confounding. vol_63d comes out with the sign of beta in OLS, in Ridge, in Lasso and in all seven windows. It is a problem with the label, not with the estimator.
• Comparing models without pairing or counting the tests. IC differences of 0.001 are noise; models are compared week by week and with the number of configurations tested in plain view (lesson 5.4.4).

Exercises

1. Add revision_analistas (the leaky feature) to the 21 and fit the Lasso with α = 1 on 2012-2018. What coefficient does it get and at what α does it enter the path? Does regularization detect that the feature is a leak?

G["revision_analistas"] = g(d["revision_analistas"])
X22 = G[ent].to_numpy()
las = Lasso(alpha=1.0).fit(X22, Y)
print(f"coefficient of the leak: {las.coef_[-1]:.1f} bp "
      f"(ret_1s: {las.coef_[0]:.1f} bp)")
al22, cam22, _ = lasso_path(X22 - X22.mean(0), Y - Y.mean(),
                            alphas=np.geomspace(80, 5, 12))
print("α at which it enters:", al22[np.argmax(cam22[-1] != 0)].round(1),
      "| features at α = 5:", int((cam22[:, -1] != 0).sum()))
G = G.drop(columns="revision_analistas")

coefficient of the leak: 69.7 bp (ret_1s: -4.0 bp)
α at which it enters: 62.2 | features at α = 5: 2

2. The panel's momentum reverses in market rebounds (momentum crash, data README). Fit OLS separately on the weeks of each regime (column regimen of panel5_coef_verdad.csv) over the whole period and compare the coefficient of mom_12_1 with the mean of the true coefficient c_mom in those weeks.

cv = pd.read_csv("panel5_coef_verdad.csv", parse_dates=["fecha"])
reg = d["fecha"].map(cv.set_index("fecha")["regimen"]).to_numpy()
cmom = d["fecha"].map(cv.set_index("fecha")["c_mom"]).to_numpy()
for r in ["calma", "estres", "rebote"]:
    m = reg == r
    b = LinearRegression().fit(G[m], y[m]).coef_[FE.index("mom_12_1")]
    print(f"{r:7s} weeks {len(np.unique(sem[m])):3d}  "
          f"OLS {b:6.2f} bp  truth {1e4 * cmom[m].mean():6.2f} bp")

calma   weeks 575  OLS   6.47 bp  truth   4.71 bp
estres  weeks  89  OLS   5.63 bp  truth   2.02 bp
rebote  weeks  65  OLS  -1.88 bp  truth -10.44 bp

3. No code: a colleague presents a Lasso on 300 automatically generated features that selects 12 and gives a mean IC of 0.08 in a five-fold k-fold. What three things would you ask for before believing it?