Learning to rank for stock selection.

Rank, do not predict

A long-short decile portfolio (lesson 5.1.3) does not use the predicted return of each stock: it uses its rank within the week. Lessons 5.5.1 to 5.5.3 train pointwise models: each row separately, with a loss that penalizes the error of each prediction, and then you sort. Learning to rank (LTR), the family of models behind search engines, trains an order directly. There are three ways to frame it:

flowchart TB
  A[One week] --> B[Pointwise]
  A --> C[Pairwise]
  A --> D[Listwise]
  B --> E[Error of each row]
  C --> F[Inverted pairs]
  D --> G[NDCG of the list]

• Pointwise: regression or classification row by row (ridge, LightGBM with squared error).
• Pairwise: a logistic loss on each pair of stocks of the same week with a different grade, which penalizes the pairs in the wrong order (RankNet; rank:pairwise in XGBoost).
• Listwise: weights each pair by how much a metric of the whole list, usually NDCG, would change if they were swapped. This is LambdaRank; with boosted trees, LambdaMART (Burges, 2010), which is what objective="lambdarank" in LightGBM and rank:ndcg in XGBoost do.

The promise is attractive: if what you trade is the order, optimizing the order should beat optimizing the squared error. The lesson tests that promise on the level's synthetic panel (README of public/curso/datos/nivel-5/), where the true alpha of each stock and week is known, and finds a trap that holds for real data.

The panel as queries

In search engine language, each query has a list of documents with an integer relevance grade. Here the query is the week, the documents are the 400 stocks and the grade, the decile of their return over the following week (0 the worst, 9 the best). The features are the panel's 21 without the leaky one, as a normal rank within each week (lesson 5.2.4); training on 2012-2018 without the last week (its label ends in 2019) and testing on 2019-2025.

flowchart TB
  A[Rows by week] --> B[group: 400 rows]
  A --> C[Grade: decile 0 to 9]
  B --> D[LGBMRanker.fit]
  C --> D
  A --> E[qid per row]
  C --> F[XGBRanker.fit]
  E --> F

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",
                    columns=["fecha", "ticker", "ret_sig_semana"])
v = pd.read_parquet("panel5_verdad.parquet",
                    columns=["fecha", "ticker", "alfa_verdadero"])
beta = pd.read_csv("panel5_activos_verdad.csv",
                   index_col="ticker")["beta"]
d = (f.merge(e, on=["fecha", "ticker"])
     .merge(v, on=["fecha", "ticker"])
     .dropna(subset=["ret_sig_semana"])
     .sort_values(["fecha", "ticker"]).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 np.nan_to_num(norm.ppf(u.to_numpy()))

X = np.column_stack([g(d[k]) for k in FE]).astype(np.float32)
y = g(d["ret_sig_semana"])                     # continuous label
r = d["ret_sig_semana"].to_numpy()
alfa = d["alfa_verdadero"].to_numpy()
b = d["ticker"].map(beta).to_numpy()
decil = (d.groupby("fecha")["ret_sig_semana"].rank(pct=True) * 10)
decil = np.minimum(decil, 9.999).astype(int).to_numpy()  # 0 to 9

fechas = np.unique(d["fecha"])
sem = np.searchsorted(fechas, d["fecha"].to_numpy())
NB = int((fechas < np.datetime64("2019-01-01")).sum())
TR = np.flatnonzero(sem < NB - 1)     # 2012-2018, without the last week
TE = np.flatnonzero(sem >= NB)        # 2019-2025
grupos = np.bincount(sem[TR])         # rows of each query
print(f"training {len(TR)} rows in {len(grupos)} weeks; "
      f"test {len(TE)} rows")
print("stocks per week:", np.unique(grupos),
      " grades:", np.bincount(decil[TR]).tolist())

training 145200 rows in 363 weeks; test 146000 rows
stocks per week: [400]  grades: [14157, 14520, 14524, 14518, 14522, 14516, 14520, 14520, 14521, 14882]

The rows have to be sorted by query: group says how many consecutive rows make up each week (363 weeks of 400) and in XGBoost qid repeats the week on each row. The grades are almost uniform, about 14,500 rows per decile (the extremes vary somewhat because of rounding ties).

NDCG: the search engine metric

For a list sorted by the model, the discounted cumulative gain of the first $k$ positions and its normalized version are

where $rel_i$ is the grade of the document in position $i$ and $IDCG@k$ the DCG of the perfect order (Järvelin and Kekäläinen, 2002). Three decisions are hidden in the formula: only the head of the list counts (the first $k$, here 40, the decile that would be bought), the position is discounted with a logarithm and the gain is exponential in the grade: decile 9 is worth 511 and decile 8, 255. Not all libraries use the same one: LightGBM uses $2^{rel}-1$ by default (label_gain), XGBoost too (ndcg_exp_gain=True, which limits grades to 31), and scikit-learn's ndcg_score uses the grade as is.

from sklearn.metrics import ndcg_score

def M(a, filas):
    """Panel vector to a weeks x 400 stocks matrix."""
    return a[filas].reshape(-1, 400)

def ndcg(p, filas, k=40, rel=None):
    """NDCG@k of each week with gain 2^rel - 1."""
    gan = 2.0 ** M(decil if rel is None else rel, filas) - 1
    orden = np.argsort(-M(p, filas), axis=1)[:, :k]
    desc = 1 / np.log2(np.arange(2, k + 2))
    dcg = (np.take_along_axis(gan, orden, 1) * desc).sum(1)
    idcg = (-np.sort(-gan, axis=1)[:, :k] * desc).sum(1)
    return dcg / idcg

azar = np.random.default_rng(0).random(len(d))
s0 = TE[:400]                                   # first week
propio = ndcg(azar, s0)[0]
sk = ndcg_score([2.0 ** decil[s0] - 1], [azar[s0]], k=40)
sk_lineal = ndcg_score([decil[s0]], [azar[s0]], k=40)
print(f"week 1: own {propio:.4f}, scikit-learn with 2^rel-1 "
      f"{sk:.4f}, scikit-learn with rel {sk_lineal:.4f}")
print(f"mean NDCG@40 of a random order (2019-2025): "
      f"{ndcg(azar, TE).mean():.4f}")

week 1: own 0.2596, scikit-learn with 2^rel-1 0.2596, scikit-learn with rel 0.5286
mean NDCG@40 of a random order (2019-2025): 0.1957

The custom function matches scikit-learn when it is given the exponential gains (0.2596 in the first week) and not when it is given the grade (0.5286): the same list has two NDCGs depending on the convention. The value of a random order, 0.1957 on average, is the reference: an NDCG only says something compared with random and with the true order.

Five models

The five are trained with the same moderate hyperparameters (200 trees, rate 0.03, large leaves) and a fixed seed; there is no search, so as not to add the selection bias of lesson 5.5.4:

• ridge on the ranks (α = 1, lesson 5.5.1) and LightGBM regression: pointwise, with the label as a normal rank.
• LightGBM LambdaRank with the default gain ($2^{rel}-1$) and with linear gain (label_gain=[0, 1, …, 9]).
• XGBoost rank:ndcg (LambdaMART) and rank:pairwise (RankNet). In XGBoost 2.0 and later, both build pairs by default with lambdarank_pair_method="topk", focused on the head of the list; exercise 3 tries random pairs.

import lightgbm as lgb
import xgboost as xgb
from sklearn.linear_model import Ridge

LGB = dict(n_estimators=200, learning_rate=0.03, num_leaves=15,
           min_child_samples=500, subsample=0.7, subsample_freq=1,
           colsample_bytree=0.7, random_state=0, n_jobs=4,
           deterministic=True, force_row_wise=True, verbose=-1)
XGB = dict(n_estimators=200, learning_rate=0.03, max_depth=4,
           min_child_weight=50, subsample=0.7, colsample_bytree=0.7,
           random_state=0, n_jobs=4, tree_method="hist")
pred = {"ridge": Ridge(alpha=1.0).fit(X[TR], y[TR]).predict(X),
        "lightgbm regresión": lgb.LGBMRegressor(**LGB)
        .fit(X[TR], y[TR]).predict(X)}
for nombre, extra in [("lambdarank 2^rel-1", {}),
                      ("lambdarank lineal",
                       {"label_gain": list(range(10))})]:
    m = lgb.LGBMRanker(objective="lambdarank", **extra, **LGB)
    pred[nombre] = m.fit(X[TR], decil[TR], group=grupos).predict(X)
for obj in ("rank:ndcg", "rank:pairwise"):
    m = xgb.XGBRanker(objective=obj, **XGB)
    pred["xgboost " + obj] = m.fit(X[TR], decil[TR],
                                   qid=sem[TR]).predict(X)
pred["oráculo"] = alfa
pred["azar"] = azar
print(list(pred))

['ridge', 'lightgbm regresión', 'lambdarank 2^rel-1', 'lambdarank lineal', 'xgboost rank:ndcg', 'xgboost rank:pairwise', 'oráculo', 'azar']

NDCG versus IC

def ic(p, filas):
    P = M(p, filas).argsort(1).argsort(1).astype(float)
    Rr = M(r, filas).argsort(1).argsort(1).astype(float)
    P -= P.mean(1, keepdims=True)
    Rr -= Rr.mean(1, keepdims=True)
    return (P * Rr).sum(1) / np.sqrt((P**2).sum(1) * (Rr**2).sum(1))

vol = X[:, FE.index("vol_63d")]
def corr_vol(p, filas):
    """Mean rank correlation with vol_63d in each week."""
    return np.mean([np.corrcoef(a, c)[0, 1] for a, c in zip(
        M(p, filas).argsort(1).argsort(1), M(vol, filas))])

pred["solo vol_63d"] = vol
tabla = pd.DataFrame({n: {"NDCG@40": ndcg(p, TE).mean(),
                          "IC": ic(p, TE).mean(),
                          "corr. con vol_63d": corr_vol(p, TE)}
                      for n, p in pred.items()}).T
tabla.round(4)

                       NDCG@40      IC  corr. con vol_63d
ridge                   0.2386  0.0180             0.0174
lightgbm regresión      0.2433  0.0125             0.0115
lambdarank 2^rel-1      0.2870  0.0086             0.8798
lambdarank lineal       0.2473  0.0101             0.2684
xgboost rank:ndcg       0.2883  0.0093             0.8748
xgboost rank:pairwise   0.2521  0.0096             0.3497
oráculo                 0.2178  0.0256            -0.3155
azar                    0.1957  0.0009             0.0001
solo vol_63d            0.2818  0.0073             0.9780

The two rankings look nothing alike. By NDCG@40, XGBoost rank:ndcg (0.2883) and LambdaRank (0.2870) win, and the oracle (oráculo), which knows the true alpha, comes almost last (0.2178), barely above random (azar). By IC, the oracle is first (0.0256), ridge second (0.0180) and the two rankers with exponential gain are the last models (0.0086 and 0.0093). The third column explains the contradiction: the scores of LambdaRank and of rank:ndcg have a rank correlation of 0.88 with vol_63d. They have learned to rank by volatility. And ranking only by vol_63d, with no model, gives an NDCG@40 of 0.2818: more than any pointwise model and more than the oracle.

Why volatility wins on NDCG

sexto = (M(vol, TE).argsort(1).argsort(1) * 6) // 400   # 0 to 5
D = M(decil, TE)
filas = {}
for q in range(6):
    m = sexto == q
    filas[q + 1] = {"en el decil alto": (D[m] == 9).mean(),
                    "en el decil bajo": (D[m] == 0).mean(),
                    "ganancia media": (2.0 ** D[m] - 1).mean()}
pd.DataFrame(filas).T.rename_axis("sexto de vol_63d").round(3)

                  en el decil alto  en el decil bajo  ganancia media
sexto de vol_63d
1                            0.052             0.057          76.680
2                            0.073             0.066          89.698
3                            0.089             0.081          97.549
4                            0.104             0.095         104.221
5                            0.127             0.121         115.219
6                            0.171             0.165         132.432

A stock in the least volatile sixth of its week ends up in the top return decile 5.2% of the time; one in the most volatile, 17.1%. But it also ends up in the bottom decile 16.5% of the time. Volatility does not predict direction; it predicts magnitude, and a stock that moves a lot reaches both extremes more often. NDCG@40 with exponential gain only looks at the head of the list and does not penalize a stock at the head ending up in decile 0 (gain 0) more than it penalizes it ending up in decile 4 (gain 15). What counts is how often the head contains deciles 8 and 9, and volatility raises that frequency: the mean gain goes from 77 in the quietest sixth to 132 in the most volatile. A model that optimizes that metric learns exactly that.

In a search engine the asymmetry makes sense: an irrelevant document at the head costs nothing. In a portfolio, a stock bought that ends up in decile 0 costs as much as one in decile 9 earns. The ranker's objective and the portfolio's are not the same.

What each model gives in a portfolio

The convention is the level's (lesson 5.1.3): each week, long the top 40 and short the bottom 40, weights of ±1/40 (gross notional of 2), simple returns and 10 bp per side on turnover. Two more columns: the portfolio beta (sum of weights times the true beta of each stock, from panel5_activos_verdad.csv; with real data it would be a beta estimated with the past) and the net Sharpe without beta, subtracting each week the portfolio beta times the average return of the 400 stocks.

def cartera(p, filas, coste_pb=10):
    """Long the top decile and short the bottom (40 and 40, weights
    ±1/40): weekly returns, turnover and portfolio beta."""
    P, Rr, B = M(p, filas), M(r, filas), M(b, filas)
    o = P.argsort(1)
    w = np.zeros_like(Rr)
    np.put_along_axis(w, o[:, -40:], 1 / 40, 1)
    np.put_along_axis(w, o[:, :40], -1 / 40, 1)
    bruto = (w * Rr).sum(1)
    rot = np.abs(np.diff(w, axis=0, prepend=0)).sum(1)
    beta_p = (w * B).sum(1)
    mercado = Rr.mean(1)
    return {"bruto": bruto, "neto": bruto - coste_pb / 1e4 * rot,
            "rot": rot, "beta": beta_p,
            "neto_sin_beta": bruto - coste_pb / 1e4 * rot
            - beta_p * mercado}

sr = lambda x: x.mean() / x.std() * np.sqrt(52)
C = {n: cartera(p, TE) for n, p in pred.items()}
eco = pd.DataFrame({n: {
    "anual bruto": c["bruto"].mean() * 52,
    "sharpe bruto": sr(c["bruto"]),
    "rotación": c["rot"].mean(), "beta": c["beta"].mean(),
    "sharpe neto": sr(c["neto"]),
    "neto sin beta": sr(c["neto_sin_beta"])} for n, c in C.items()}).T
eco.round(3)

                       anual bruto  sharpe bruto  rotación   beta  sharpe neto  neto sin beta
ridge                        0.128         1.364     2.286 -0.027        0.096          0.188
lightgbm regresión           0.101         1.265     2.622  0.000       -0.434         -0.266
lambdarank 2^rel-1           0.086         0.618     1.297  0.470        0.136         -0.542
lambdarank lineal            0.127         1.476     2.575  0.165       -0.077         -0.276
xgboost rank:ndcg            0.099         0.744     1.272  0.473        0.247         -0.386
xgboost rank:pairwise        0.113         1.309     2.552  0.133       -0.234         -0.391
oráculo                      0.167         1.676     2.648 -0.218        0.295          0.745
azar                        -0.006        -0.110     3.599 -0.005       -3.411         -3.410
solo vol_63d                 0.104         0.654     0.438  0.617        0.510          0.034

The figure sorts the three Sharpe columns by model. Three readings:

• Gross, LambdaRank with linear gain (1.48) and ridge (1.36) are at the top, close to the oracle (1.68). The exponential rankers give half (0.62 and 0.74).
• Net, the exponential rankers seem to fare better (0.14 and 0.25 versus 0.10 for ridge): they turn over half as much (1.30 and 1.27 weekly versus 2.29), because volatility changes slowly. Ranking only by vol_63d turns over 44% and gives the best net of all (0.51).
• Without beta, the mirage falls apart. The portfolios of the exponential rankers have a beta of 0.47 (0.62 for volatility alone): they are long volatile stocks and short quiet ones, a bet on the market, which rises in 2019-2025. Removing beta, LambdaRank gives −0.54 and rank:ndcg, −0.39; volatility alone, 0.03. Ridge, with beta −0.03, keeps 0.19, and the oracle 0.75 (its beta is negative, −0.22, because of the planted low-volatility premium).

The header figure sums up the lesson: on this panel, the more NDCG@40, the worse the economic result.

Does any of them beat ridge?

The fair comparison is paired: in the same week, the difference in beta-free net return between each model and ridge, with its standard error.

base = C["ridge"]["neto_sin_beta"]
for n in ["lightgbm regresión", "lambdarank 2^rel-1",
          "lambdarank lineal", "xgboost rank:ndcg",
          "xgboost rank:pairwise"]:
    dif = C[n]["neto_sin_beta"] - base
    t = dif.mean() / dif.std() * np.sqrt(len(dif))
    print(f"{n:22s} {dif.mean() * 52:+.2%} annual  t = {t:+.2f}")

lightgbm regresión     -3.69% annual  t = -1.29
lambdarank 2^rel-1     -6.57% annual  t = -1.45
lambdarank lineal      -3.77% annual  t = -1.20
xgboost rank:ndcg      -5.07% annual  t = -1.16
xgboost rank:pairwise  -4.67% annual  t = -1.55

They all lose to ridge, by −3.7% to −6.6% annualized, with t between −1.2 and −1.6: no difference is significant on its own, but they all go in the same direction, and ridge is the cheapest and the easiest to explain. That the best one is the baseline does not make it a strategy: with six models tried over the same period, the best one is also selected.

def psr(s, umbral, t, asim, curt):
    den = np.sqrt(1 - asim * s + (curt - 1) / 4 * s**2)
    return norm.cdf((s - umbral) * np.sqrt(t - 1) / den)

def max_esperado(n):
    gm = np.euler_gamma
    return ((1 - gm) * norm.ppf(1 - 1 / n)
            + gm * norm.ppf(1 - 1 / (n * np.e)))

modelos = ["ridge", "lightgbm regresión", "lambdarank 2^rel-1",
           "lambdarank lineal", "xgboost rank:ndcg",
           "xgboost rank:pairwise"]
S = pd.DataFrame({n: C[n]["neto_sin_beta"] for n in modelos})
srs = S.mean() / S.std()
k = srs.idxmax()
T, asim, curt = len(S), S[k].skew(), S[k].kurt() + 3
sr0 = np.sqrt(srs.var()) * max_esperado(len(modelos))
print(f"best: {k}, annualized Sharpe {srs[k] * np.sqrt(52):.2f}")
print(f"PSR(0) {psr(srs[k], 0, T, asim, curt):.3f}; "
      f"DSR with N = {len(modelos)}: "
      f"{psr(srs[k], sr0, T, asim, curt):.3f}")

best: ridge, annualized Sharpe 0.19
PSR(0) 0.692; DSR with N = 6: 0.359

Ridge's Sharpe of 0.19 has a PSR(0) of 0.692 and a Deflated Sharpe of 0.359 with N = 6: there is no evidence that any of these portfolios has a positive return after costs, beta and selection (lessons 2.5.4 and 5.4.4). It is what the panel's README anticipates: with 10 bp per side, not even the oracle gets past a modest net Sharpe, and the weekly turnover takes almost everything (lesson 5.1.3 smooths the signal to turn over less).

When learning to rank makes sense

None of the above says that LTR is useless in finance; it says that a search engine's objective is not a portfolio's. Poh, Lim, Zohren and Roberts (2021) incorporate pairwise and listwise rankers into a cross-sectional momentum strategy and, on their data, report Sharpes about three times higher than those of ranking with heuristics or with the outputs of regressions and classifiers. It is a result with their data and their validation; when replicating it, the first thing is to check this lesson's trap: how much of the improvement is magnitude and beta. For a ranker to serve a long-short portfolio it helps to have:

• Gains that weigh both extremes equally. Linear gain, symmetric grades or two rankers, one for each side; exercise 1 shows that training a second LambdaRank for the short side is not enough if the gain is still exponential.
• Labels without the magnitude you do not want. Returns residual to beta and sector, or standardized by a volatility known at the date, so that each stock's dispersion does not decide the grade.
• Portfolio metrics. Select the model by IC, net Sharpe and beta, never by NDCG alone, and compare with the pointwise baseline over the same weeks.
• A horizon and a head size that match the portfolio. lambdarank_truncation_level (30 by default in LightGBM) and the $k$ of NDCG should be the number of stocks bought.

Common mistakes

• Reading NDCG without the random and oracle references. 0.287 looks like much more than 0.218, and it is the difference between betting on volatility and knowing the alpha.
• Mixing up NDCG conventions. The same list gives 0.26 or 0.53 depending on the gain; always compare with the same function.
• Rows not sorted by query. group and qid assume contiguous rows for each week; if they are not, the ranker compares stocks from different weeks without warning.
• Grades that mix magnitude and direction. Deciles of raw return reward volatility under an asymmetric metric.
• Attributing to selection what is beta. A long-short of volatile against quiet stocks wins when the market rises.
• Not comparing with the pointwise baseline. Ridge beats all the rankers on this panel.

Exercises

1. Train a second LambdaRank with inverted grades (9 − decil) to choose the stocks to short, and compare the correlation of its scores with vol_63d and its IC with those of the LambdaRank for the long side.

m = lgb.LGBMRanker(objective="lambdarank", **LGB)
corto = -m.fit(X[TR], 9 - decil[TR], group=grupos).predict(X)
largo = pred["lambdarank 2^rel-1"]
for n, p in [("long", largo), ("short", corto)]:
    print(f"{n}: corr. with vol_63d {corr_vol(p, TE):+.2f}, "
          f"IC {ic(p, TE).mean():+.4f}")

long: corr. with vol_63d +0.88, IC +0.0086
short: corr. with vol_63d -0.89, IC -0.0058

2. Compute NDCG@40 with linear gain and the NDCG of the whole list (k = 400) with ndcg_score for ridge, the two LambdaRanks, the oracle, vol_63d alone and random. Which version separates the oracle from random best?

def ndcg_sk(p, k=None):
    return np.mean([ndcg_score([dd], [pp], k=k) for dd, pp
                    in zip(M(decil, TE), M(p, TE))])

lineal = {}
for n in ["ridge", "lambdarank 2^rel-1", "lambdarank lineal",
          "oráculo", "solo vol_63d", "azar"]:
    lineal[n] = {"NDCG@40 2^rel-1": ndcg(pred[n], TE).mean(),
                 "NDCG@40 lineal": ndcg_sk(pred[n], 40),
                 "NDCG@400 lineal": ndcg_sk(pred[n])}
pd.DataFrame(lineal).T.round(4)

                    NDCG@40 2^rel-1  NDCG@40 lineal  NDCG@400 lineal
ridge                        0.2386          0.5139           0.8733
lambdarank 2^rel-1           0.2870          0.5023           0.8712
lambdarank lineal            0.2473          0.5147           0.8733
oráculo                      0.2178          0.5210           0.8748
solo vol_63d                 0.2818          0.5026           0.8709
azar                         0.1957          0.5007           0.8698

3. Train XGBoost's rank:pairwise with random pairs over the whole list (lambdarank_pair_method="mean", one pair per stock) instead of the default topk. Does the volatility bias disappear? Does the portfolio improve?

m = xgb.XGBRanker(objective="rank:pairwise",
                  lambdarank_pair_method="mean",
                  lambdarank_num_pair_per_sample=1, **XGB)
p = m.fit(X[TR], decil[TR], qid=sem[TR]).predict(X)
c = cartera(p, TE)
print(f"NDCG@40 {ndcg(p, TE).mean():.4f}, IC {ic(p, TE).mean():.4f}, "
      f"corr. with vol_63d {corr_vol(p, TE):+.3f}")
print(f"gross Sharpe {sr(c['bruto']):.2f}, net {sr(c['neto']):.2f}, "
      f"turnover {c['rot'].mean():.2f}")

NDCG@40 0.2191, IC 0.0036, corr. with vol_63d +0.027
gross Sharpe 0.48, net -1.88, turnover 2.91