Data drift and concept drift: detecting them and what to do.

Two drifts, and only one hurts

A supervised model learns the relationship between some features $X$ and a target $y$ from a sample. What it learns is an approximation to $P(y \mid X)$, and what it uses to estimate it are examples drawn from $P(X)$. Either of the two can change once the model is already in production, and the literature separates them (Gama et al., 2014; Quiñonero-Candela et al., 2009):

Concept drift, in turn, can be abrupt (a regulatory change, a competitor arbitraging the signal), gradual (an anomaly that is published and slowly fades) or recurrent (calm and stress regimes that come back).

flowchart TB
  D[Drift] --> X["Data: P(X)"]
  D --> C["Concept: P(y|X)"]
  C --> B[Abrupt]
  C --> G[Gradual]
  C --> R[Recurrent]

Data drift is the one you can see: it is enough to compare today's feature distribution with the training one, and that can be done every day, without waiting for any label. Concept drift is the one that hurts, and it cannot be seen by looking at the features: it only shows up in the results, once the labels the model predicted wrongly have arrived. In real data you never know which of the two has happened or when. The level's deriva5_* data do tell you.

The lab: six segments with known dates

deriva5_muestras.parquet has 30 fictitious assets over 3,018 NYSE sessions (2014-2025), six standardized features and the next day's return already aligned. The target is linear in the features with the signal-to-noise ratio of daily returns: a true IC of 0.10, an $R^2$ of 1%. The level's README publishes the segments:

The protagonist is a fixed model: an OLS estimated on 2014-2016 and never retrained. Its quality is measured as throughout the level, with the daily IC: the cross-sectional correlation between its prediction and the next day's return of the 30 assets (lesson 5.1.3).

import numpy as np
import pandas as pd

d = pd.read_parquet("deriva5_muestras.parquet")
seg = pd.read_csv("deriva5_verdad.csv", parse_dates=["desde"])
F = [f"f{k}" for k in range(1, 7)]
fechas = pd.DatetimeIndex(np.sort(d["fecha"].unique()))
n, NA = len(fechas), d["activo"].nunique()
X = d[F].to_numpy(float).reshape(n, NA, 6)  # session × asset × f
Y = d["objetivo"].to_numpy(float).reshape(n, NA)
s = np.searchsorted(seg["desde"], fechas, side="right") - 1

base = s == 0                                # 2014-2016
X1 = np.column_stack([np.ones(base.sum() * NA),
                      X[base].reshape(-1, 6)])
beta = np.linalg.lstsq(X1, Y[base].ravel(), rcond=None)[0]
print("fixed model coefficients (bp):",
      np.round(beta[1:] * 1e4, 1))

def ic_diario(P, Y):
    """Cross-sectional correlation between P and Y in each session."""
    P = P - P.mean(1, keepdims=True)
    Yc = Y - Y.mean(1, keepdims=True)
    return ((P * Yc).sum(1)
            / np.sqrt((P**2).sum(1) * (Yc**2).sum(1)))

ic = pd.Series(ic_diario(X @ beta[1:], Y), index=fechas)
tabla = seg[["segmento", "desde", "tipo"]].copy()
tabla["desde"] = tabla["desde"].dt.date
tabla["ic_fijo"] = ic.groupby(s).mean().round(3).to_numpy()
tabla

fixed model coefficients (bp): [ 8.9  6.1 -4.5 -0.  -0.   1. ]
   segmento       desde               tipo  ic_fijo
0         0  2014-01-02               base    0.120
1         1  2017-01-03              datos    0.117
2         2  2019-07-01           concepto   -0.019
3         3  2021-01-04  datos_irrelevante   -0.026
4         4  2022-07-01   concepto_gradual   -0.035
5         5  2023-07-03            estable   -0.051

The OLS recovers the base relationship with its sampling error (8.9, 6.1 and −4.5 bp against 8, 5 and −4) and zeros for the features with no effect. Its mean IC does not move with the 2017 data drift (0.117 against 0.120) and sinks with the 2019 concept drift (−0.019). From then on it does not recover: in 2021 the drift of f4 is irrelevant to it, and the gradual drift of 2022 also takes away the f2 part (−0.051 at the end).

Data drift: PSI and KS

The Population Stability Index (PSI) comes from credit risk. The reference data are cut into $B$ bins (here, deciles), the fraction $p_i$ of the reference and $q_i$ of the current period falling in each one is measured, and the sum is taken:

It is a symmetrized Kullback-Leibler divergence over the histogram: zero if the two histograms coincide, and it grows with any difference, in mean, dispersion or shape. The usual rule of thumb reads a PSI below 0.1 as stable and above 0.25 as strong drift.

def psi(ref, cur, nb=10):
    """Population Stability Index with deciles of the reference."""
    cortes = np.quantile(ref, np.linspace(0, 1, nb + 1)[1:-1])
    p = np.bincount(np.searchsorted(cortes, ref), minlength=nb)
    q = np.bincount(np.searchsorted(cortes, cur), minlength=nb)
    p = np.maximum(p / len(ref), 1e-6)
    q = np.maximum(q / len(cur), 1e-6)
    return float(np.sum((q - p) * np.log(q / p)))

ref = X[base]
filas = {}
for k in range(1, 6):
    cur = X[s == k]
    filas[f"{k} {seg['tipo'][k]}"] = [
        psi(ref[..., j].ravel(), cur[..., j].ravel())
        for j in range(6)]
pd.DataFrame(filas, index=F).T.round(2)

                       f1   f2   f3    f4   f5   f6
1 datos              0.97  0.0  0.5  0.00  0.0  0.0
2 concepto           0.00  0.0  0.0  0.00  0.0  0.0
3 datos_irrelevante  0.00  0.0  0.0  1.23  0.0  0.0
4 concepto_gradual   0.00  0.0  0.0  1.27  0.0  0.0
5 estable            0.00  0.0  0.0  1.28  0.0  0.0

The PSI sees exactly the planted data drifts: 0.97 in f1 and 0.50 in f3 in segment 1 (the doubled dispersion of f3 counts as well), 1.23 in f4 from 2021, and zeros everywhere else. Above all, zeros in segment 2, the concept drift one, and in the gradual one. The figure at the top puts the two things together: the PSI jumps in 2017 and in 2021 without the IC noticing, and the IC sinks in July 2019 without any PSI moving.

The classic alternative is the two-sample Kolmogorov-Smirnov test (scipy.stats.ks_2samp), which compares the empirical distribution functions and gives a p-value. Before using either of them as an alarm you need to know what they give when nothing is happening: each month of 2016 is compared with 2014-2015, with no planted change.

from scipy import stats

ref = X[fechas.year <= 2015].reshape(-1, 6)   # 2014-2015
meses = pd.period_range("2016-01", "2016-12", freq="M")
psis, pvals = [], []
for m in meses:
    cur = X[fechas.to_period("M") == m].reshape(-1, 6)
    for j in range(6):
        psis.append(psi(ref[:, j], cur[:, j]))
        pvals.append(stats.ks_2samp(ref[:, j], cur[:, j]).pvalue)
psis, pvals = np.array(psis), np.array(pvals)
print(f"72 tests (12 months × 6 features) with no drift")
print(f"PSI: median {np.median(psis):.3f}, maximum {psis.max():.3f}")
print(f"KS with p < 0.05: {np.mean(pvals < 0.05):.0%} "
      f"(expected: 5%)")
ks_ano = [stats.ks_2samp(ref[:, j],
                         X[fechas.year == 2016][..., j].ravel())
          for j in range(6)]
print("all of 2016, KS p-value:",
      [f"{r.pvalue:.3f}" for r in ks_ano])
print("D statistic:", [round(float(r.statistic), 3)
                       for r in ks_ano])

72 tests (12 months × 6 features) with no drift
PSI: median 0.018, maximum 0.055
KS with p < 0.05: 33% (expected: 5%)
all of 2016, KS p-value: ['0.000', '0.146', '0.009', '0.000', '0.315', '0.018']
D statistic: [0.029, 0.016, 0.023, 0.04, 0.013, 0.022]

Without drift, the monthly PSI stays at noise levels (median 0.018, maximum 0.055), well below 0.25. KS does not: 33% of the 72 tests come out with $p < 0.05$, six times what is expected, and with all of 2016 against 2014-2015 four of the six features “change” with $p < 0.05$ (three with $p \le 0.009$). There are two reasons. The features are AR(1) with $\rho = 0.5$ within each asset: the rows are not independent, the effective sample size is several times smaller than the number of rows and the p-value comes out too small. And with thousands of rows, tiny differences are significant: the D statistics are between 0.013 and 0.040, differences of a few hundredths in the distribution function.

Multivariate drift: adversarial validation

PSI and KS look at each feature separately. A drift in the relationship between features (two that stop being correlated) cannot be seen that way. Adversarial validation looks at all of them at once: a classifier is trained to distinguish the reference data from the current data. If it cannot (AUC ≈ 0.5), the two samples are indistinguishable; if it can, the classifier's importance says which features have changed. Rabanser et al. (2019) compare several drift detectors on several datasets: the best in their study is a two-sample test on a reduced representation, and classifiers that separate domains turn out to be useful for describing the drift and seeing whether it is harmful.

Two details decide whether the AUC means anything. The folds must be whole months (rows from neighboring days are almost identical, as in lesson 5.4.1) and the AUC is computed in each fold and averaged. If the predictions of all the folds are pooled and a single AUC is computed, it comes out below 0.5 with no drift: each fold has a different class proportion and the classifier, with nothing to learn, predicts that proportion, which is the wrong way round in the test fold.

import lightgbm as lgb
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import StratifiedGroupKFold

mref = fechas.year == 2016                 # the last year seen
mes = np.repeat((fechas.year * 100 + fechas.month).to_numpy(), NA)
mes = mes.reshape(n, NA)

def adversaria(m):
    """AUC of a classifier that separates 2016 from sessions m."""
    A = np.vstack([X[mref].reshape(-1, 6), X[m].reshape(-1, 6)])
    c = np.r_[np.zeros(mref.sum() * NA), np.ones(m.sum() * NA)]
    g = np.r_[mes[mref].ravel(), mes[m].ravel()]
    aucs, imp = [], np.zeros(6)
    cv = StratifiedGroupKFold(4, shuffle=True, random_state=0)
    for tr, te in cv.split(A, c, g):          # whole months
        clf = lgb.LGBMClassifier(
            n_estimators=60, num_leaves=8, min_child_samples=100,
            random_state=0, n_jobs=4, deterministic=True,
            force_row_wise=True, verbose=-1).fit(A[tr], c[tr])
        p = clf.predict_proba(A[te])[:, 1]
        aucs.append(roc_auc_score(c[te], p))  # AUC of each fold
        imp += clf.booster_.feature_importance("gain")
    return np.mean(aucs), F[int(imp.argmax())], imp.max() / imp.sum()

semestre = fechas.year.astype(str) + np.where(
    fechas.month <= 6, "-S1", "-S2")
adv = {}
for sm in semestre.unique():
    if sm.startswith("2016"):
        continue
    auc, top, peso = adversaria(semestre == sm)
    adv[sm] = {"auc": auc, "feature": top, "weight": peso,
               "ic_fijo": ic[semestre == sm].mean()}
adv = pd.DataFrame(adv).T.astype({"auc": float, "weight": float,
                                  "ic_fijo": float})
adv.loc[["2014-S1", "2015-S2", "2017-S1", "2019-S1", "2019-S2",
         "2020-S2", "2021-S1", "2023-S1", "2025-S2"]].round(3)

           auc feature  weight  ic_fijo
2014-S1  0.493      f3   0.203    0.138
2015-S2  0.518      f4   0.198    0.103
2017-S1  0.826      f1   0.535    0.142
2019-S1  0.824      f1   0.533    0.123
2019-S2  0.519      f2   0.198   -0.031
2020-S2  0.496      f1   0.185   -0.025
2021-S1  0.782      f4   0.905   -0.029
2023-S1  0.780      f4   0.916   -0.026
2025-S2  0.790      f4   0.919   -0.027

In the half-years without drift the AUC stays between 0.49 and 0.52, what two identical samples give. In 2017-2019 it rises to 0.83, with f1 as the most important feature (53% of the gain). From 2021 it stays at 0.78-0.79, with f4 (92%). And in 2019-S2 and 2020-S2, after the concept change, it goes back to 0.52 and 0.50 while the IC of the fixed model is negative (−0.031 and −0.025). Adversarial validation answers its own question very well, whether today's data look like the training data, and that question is not the one that matters for the P&L.

So what is monitoring the features good for? Diagnosing (when performance drops, knowing whether the features changed helps to understand why), warning about data errors (a vendor that changes a unit, a column that arrives empty) and monitoring while the labels have not arrived. With a 21-day target, today's performance is known within a month; the features, today.

Monitoring performance: moving average, CUSUM and Page-Hinkley

Concept drift can only be seen in performance. The daily IC is a noisy series: with 30 assets and a true IC of 0.12, its standard deviation is 0.18 (almost $1/\sqrt{29}$), so a day with a negative IC says nothing. The simplest detector, the one in the level's README, is the 63-session moving average of the IC with an alarm when it falls below zero.

cambio = pd.Timestamp("2019-07-01")
mm = ic.rolling(63).mean()
fecha_alarma = mm[(mm.index >= cambio) & (mm < 0)].index[0]
retraso = ((fechas >= cambio) & (fechas <= fecha_alarma)).sum()
print("first moving average < 0:", fecha_alarma.date(),
      f"({retraso} sessions after the change)")
antes = mm["2017-01-01":"2019-06-28"]
print(f"minimum before the change: {antes.min():.3f}")

first moving average < 0: 2019-08-29 (43 sessions after the change)
minimum before the change: 0.071

It alarms on 2019-08-29, 43 sessions after the change, with no false alarms before (the moving average did not fall below 0.071 in 2017-2019). It is easy to explain and slow: a 63-day window needs half of its content to be after the change for the mean to cross zero.

Page's (1954) CUSUM is the classic sequential detector. For a drop in the mean from $\mu_0$:

$S_t$ accumulates how far the IC stays below $\mu_0 - k$ and resets to zero when the model is doing well. The slack $k$ is usually set at half of the jump to be detected (here, from 0.12 to 0: $k = 0.06$) and the threshold $h$ decides the balance between false alarms and delay. It is the same recursion as the CUSUM filter of lesson 5.3.5, which there sampled prices and here monitors a model. Page-Hinkley is its variant without $\mu_0$: it compares each data point with the mean since the last reset.

x = ic.to_numpy()
t0 = np.searchsorted(fechas, pd.Timestamp("2017-01-03"))
mu0 = ic[base].mean()

def cusum(x, mu0, k, h):
    """One-sided Page CUSUM for a drop in the mean."""
    S, alarmas = 0.0, []
    for t, v in enumerate(x):
        S = max(0.0, S + (mu0 - k) - v)
        if S > h:
            alarmas.append(t)
            S = 0.0                       # reset after the alarm
    return alarmas

def page_hinkley(x, delta, lam):
    """Page-Hinkley: drop against the mean since the reset."""
    m = M = suma = 0.0
    cnt, alarmas = 0, []
    for t, v in enumerate(x):
        cnt += 1
        suma += v
        m += suma / cnt - v - delta
        M = min(M, m)
        if m - M > lam:
            alarmas.append(t)
            m = M = suma = 0.0
            cnt = 0
    return alarmas

def resumen(alarmas):
    f = fechas[t0:][alarmas]
    falsas = int((f < cambio).sum())
    tras = f[f >= cambio]
    r = int(((fechas >= cambio) & (fechas <= tras[0])).sum())
    return falsas, tras[0].date(), r

print(f"mean reference IC (2014-2016): {mu0:.3f}")
for h in [0.5, 1.0, 2.0, 3.0]:
    fa, f, r = resumen(cusum(x[t0:], mu0, 0.06, h))
    print(f"CUSUM h={h}: {fa:2d} false alarms 2017-2019, "
          f"first after the change {f} ({r} sessions)")
for lam in [1.0, 2.0, 3.0]:
    fa, f, r = resumen(page_hinkley(x[t0:], 0.03, lam))
    print(f"Page-Hinkley lambda={lam}: {fa:2d} false alarms, "
          f"{f} ({r} sessions)")

mean reference IC (2014-2016): 0.120
CUSUM h=0.5: 19 false alarms 2017-2019, first after the change 2019-07-08 (5 sessions)
CUSUM h=1.0:  3 false alarms 2017-2019, first after the change 2019-07-16 (11 sessions)
CUSUM h=2.0:  0 false alarms 2017-2019, first after the change 2019-07-19 (14 sessions)
CUSUM h=3.0:  0 false alarms 2017-2019, first after the change 2019-08-01 (23 sessions)
Page-Hinkley lambda=1.0:  8 false alarms, 2019-07-11 (8 sessions)
Page-Hinkley lambda=2.0:  1 false alarms, 2019-07-18 (13 sessions)
Page-Hinkley lambda=3.0:  0 false alarms, 2019-07-29 (20 sessions)

The CUSUM with $h = 2$ alarms on 2019-07-19, after 14 sessions, with no false alarm in the previous two and a half years, which also contain the 2017 data drift. With $h = 0.5$ it alarms after 5 sessions, at the cost of 19 false alarms over that period, about seven per year. Page-Hinkley behaves similarly: with $\lambda = 3$, 20 sessions with no false alarms. The top half of the figure follows the CUSUM statistic (without reset) in 2019: flat until July and, from the change on, an almost straight rise.

False alarms against delay

A detector is not judged by one episode. The 2019 change is a single realization, and 14 sessions can be luck. To measure the trade-off, many histories are simulated: 2,000 ten-year paths with ICs drawn at random from the days without change (2017-2019-S1) to count false alarms, and one-year paths with ICs from after the change (2019-S2 to 2020) to measure the delay.

tc = np.searchsorted(fechas, cambio)
t21 = np.searchsorted(fechas, pd.Timestamp("2021-01-01"))
sano, roto = x[t0:tc], x[tc:t21]     # IC before and after
rng = np.random.default_rng(0)
NS, T = 2000, 2520                   # 2,000 paths of 10 years
a = rng.choice(sano, (NS, T))
b = rng.choice(roto, (NS, 252))

def sim_cusum(h, k=0.06):
    S, fa = np.zeros(NS), np.zeros(NS)
    for t in range(T):
        S = np.maximum(0, S + (mu0 - k) - a[:, t])
        fa += S > h
        S[S > h] = 0
    S, r = np.zeros(NS), np.full(NS, 252)
    for t in range(252):
        S = np.maximum(0, S + (mu0 - k) - b[:, t])
        r[(S > h) & (r == 252)] = t + 1
    return fa.mean() / 10, np.median(r)

def sim_media(w):
    c = np.cumsum(a, 1)
    bajo = (c[:, w:] - c[:, :-w]) / w < 0
    fa = (bajo[:, 1:] & ~bajo[:, :-1]).sum(1).mean() / 10
    c = np.cumsum(np.hstack([np.zeros((NS, 1)), a[:, -w:], b]), 1)
    tras = ((c[:, w:] - c[:, :-w]) / w)[:, 1:] < 0
    r = np.where(tras.any(1), tras.argmax(1) + 1, 252)
    return fa, np.median(r)

filas = [("CUSUM", h, *sim_cusum(h)) for h in [1, 1.5, 2, 2.5, 3]]
filas += [("moving average", w, *sim_media(w)) for w in [21, 42, 63]]
pd.DataFrame(filas, columns=["detector", "parameter",
                             "false/year", "median delay"]
             ).round(3)

         detector  parameter  false/year  median delay
0           CUSUM        1.0       1.030          11.0
1           CUSUM        1.5       0.193          17.0
2           CUSUM        2.0       0.038          23.0
3           CUSUM        2.5       0.007          30.0
4           CUSUM        3.0       0.001          36.0
5  moving average       21.0       0.204          18.0
6  moving average       42.0       0.002          36.0
7  moving average       63.0       0.000          54.0

The typical delay of the CUSUM with $h = 2$ is 23 sessions, not 14, with 0.04 false alarms per year (one every 25 years). The 63-session moving average takes a median of 54: the 43 in the data were also somewhat lucky. The figure puts the two families on the same plane. At the same false alarm rate, the CUSUM arrives a few sessions earlier: with 0.2 per year, 17 sessions against 18; with about 0.002-0.007, 30 against 36. The advantage exists (the CUSUM is optimal for this problem in Lorden's sense), but it is modest. The important decision is where to sit on the curve: how much a false alarm costs (switching off a model that works) against how much each day of delay costs (trading with a broken model).

Regimes: the two-state filter

Lesson 2.1.7 detected regimes with a hidden Markov model and distinguished between filtered probabilities (with the information up to each day) and smoothed ones (with the whole sample, including the future). The same idea monitors a model: two states, “works” (mean IC $\mu_0 = 0.12$) and “does not work” (mean IC 0), with the same deviation, a probability of switching state of 0.002 per day and the filter recursion:

where $A$ is the transition matrix and $\phi(x_t)$ the density of today's IC in each state. The parameters are fixed in advance, with 2014-2016, instead of estimating them on the whole sample: estimating them with data from 2019-2025 would be looking at the future even though the filter is used.

from scipy.stats import norm

medias = np.array([mu0, 0.0])     # the model works / does not work
sd = ic[base].std()
p = 0.998                         # expected duration: 500 sessions
A = np.array([[p, 1 - p], [1 - p, p]])
pi, prob = np.array([1.0, 0.0]), []
for v in x[t0:]:
    pi = (pi @ A) * norm.pdf(v, medias, sd)  # predict and update
    pi = pi / pi.sum()
    prob.append(pi[1])
prob = pd.Series(prob, index=fechas[t0:])
print(f"std of the daily IC: {sd:.3f}")
print("maximum P(broken) before the change:",
      round(prob[:"2019-06-28"].max(), 2))
print("first P(broken) > 0.9:",
      prob[(prob.index >= cambio) & (prob > 0.9)].index[0].date())
print(prob.resample("QE").mean()["2019":"2020"].round(2)
      .to_string())

std of the daily IC: 0.184
maximum P(broken) before the change: 0.63
first P(broken) > 0.9: 2019-07-19
2019-03-31    0.02
2019-06-30    0.01
2019-09-30    0.83
2019-12-31    0.98
2020-03-31    0.98
2020-06-30    0.98
2020-09-30    0.97
2020-12-31    0.99
Freq: QE-DEC

The probability that the model is broken exceeds 0.9 on 2019-07-19, the same day as the CUSUM. It is no coincidence: with two Gaussian states with known means, the log-likelihood ratio of each day is linear in $x_t$ with a slope proportional to the difference of means, and the CUSUM with $k$ at the midpoint accumulates exactly that ratio. The filter adds the prior probability of change and gives a figure that can be read directly. Its closest false alarm before the change reached 0.63, the bottom half of the figure.

Recurrent regimes call for a different response. In the level's panel momentum reverses in every rebound (lesson 5.1.1): a detector that alarms once the rebound has already started arrives late to a regime that lasts a few months, and retraining with rebound data teaches the model a relationship that is going to reverse again. If the regimes come back, it is better for the model to see them as an observable feature and learn the interaction (lesson 5.10.1); drift detection is for what does not come back.

What to do: retrain, reweight, switch off

An alarm does not say what to do. First you diagnose: if the features also changed (PSI, adversarial validation), it may be data drift with a real effect on the model, a data error or both; if not, it is concept drift.

flowchart TB
  M[IC monitor] --> A{Alarm?}
  A -- no --> M
  A -- yes --> D{High PSI or AUC?}
  D -- no --> C[Concept]
  D -- yes --> E[Data or errors]
  C --> R[Switch off, retrain]
  E --> R

The possible responses, compared in a walk-forward that retrains on the first day of each month with the targets already known:

• Retrain with an expanding window: the whole past.
• Retrain with a rolling window: only the last year (252 sessions) or half year (126).
• Reweight: the whole past with weights that halve every 126 sessions.
• The oracle: the true coefficients of each day, the ceiling.

def ols(i0, i1, semivida=None):
    """OLS with sessions i0..i1-1 (exponential weights)."""
    A = np.column_stack([np.ones((i1 - i0) * NA),
                         X[i0:i1].reshape(-1, 6)])
    y = Y[i0:i1].ravel()
    if semivida:
        edad = np.repeat(np.arange(i1 - i0)[::-1], NA)
        w = np.sqrt(0.5 ** (edad / semivida))
        A, y = A * w[:, None], y * w
    return np.linalg.lstsq(A, y, rcond=None)[0][1:]

mes_n = (fechas.year * 12 + fechas.month).to_numpy()
reentreno = {t for t in range(t0, n) if mes_n[t] != mes_n[t - 1]}

def walk_forward(entrenar):
    """Retrains on the first day of each month with what is known."""
    P, b = np.zeros((n, NA)), None
    for t in range(t0, n):
        if t in reentreno:
            b = entrenar(t - 1)      # targets known up to t
        P[t] = X[t] @ b
    return P

b_verdad = pd.read_csv("deriva5_coef_verdad.csv").iloc[:, 1:]
esquemas = {
    "fixed (2014-2016)": lambda t: beta[1:],
    "expanding": lambda t: ols(0, t),
    "rolling 252": lambda t: ols(t - 252, t),
    "rolling 126": lambda t: ols(t - 126, t),
    "weights, half-life 126": lambda t: ols(0, t, 126),
}
pred = {k: walk_forward(f) for k, f in esquemas.items()}
pred["oracle (true b)"] = np.einsum(
    "tik,tk->ti", X, b_verdad.to_numpy())
res = {}
for k, P in pred.items():
    v = pd.Series(ic_diario(P[t0:], Y[t0:]), index=fechas[t0:])
    res[k] = {"2017-2025": v.mean(),
              "2019-S2 to 2020": v["2019-07":"2020-12"].mean(),
              "2022-S2 to 2023-S1": v["2022-07":"2023-06"].mean()}
pd.DataFrame(res).T.round(3)

                        2017-2025  2019-S2 to 2020  2022-S2 to 2023-S1
fixed (2014-2016)           0.007           -0.019              -0.035
expanding                   0.052            0.001               0.027
rolling 252                 0.101            0.055               0.105
rolling 126                 0.101            0.066               0.098
weights, half-life 126      0.101            0.054               0.104
oracle (true b)             0.111            0.091               0.108

Over the whole 2017-2025 period the fixed model barely adds anything (IC 0.007), the expanding window reaches 0.052 and the three variants that forget the past, 0.101, close to the oracle (0.111). The difference lies in the 18 months that follow the concept change: the expanding window gives 0.001, because up to 2019 it accumulates five years of a positive f1 effect; the 126-session rolling window gives 0.066 and the 252-session one, 0.055. In the gradual period of 2022-2023 all the adaptive schemes follow the oracle (0.098-0.105 against 0.108): a change spread over twelve months is easy to follow with a one-year window.

This is the opposite of what lesson 5.4.3 found on the weekly panel, where the expanding window won, and there is no contradiction. There the effects are small and oscillate without a trend; here the change is large relative to the noise and permanent. The optimal window depends on that (Pesaran and Timmermann, 2007, cited in lesson 5.4.3).

The figure follows the coefficient of f1 estimated each month. The 126-session rolling window crosses zero at the end of 2019 and the 252-session one in January 2020; afterwards they oscillate around −8 with their estimation error. The expanding window was already falling before the change (from 9 to 6.5 in 2017-2019: the shift of f1 in 2017, with no change in the target's mean, attenuates the slope of an OLS with a single constant for both periods) and does not turn negative until November 2024.

The combination with the detector: the fixed model until the CUSUM alarm; after that, switched off for a few sessions (no position) and retrained only with the data after the alarm.

alarma = t0 + cusum(x[t0:], mu0, 0.06, 2.0)[0]
print("CUSUM alarm (h = 2):", fechas[alarma].date())

def tras_alarma(espera):
    """Fixed model until the alarm; switched off for `espera`
    sessions; then OLS only with data since the alarm (monthly)."""
    P, b = np.zeros((n, NA)), beta[1:]
    for t in range(t0, n):
        if alarma <= t < alarma + espera:
            b = np.zeros(6)
        elif t >= alarma + espera and (t in reentreno
                                       or t == alarma + espera):
            b = ols(alarma, t - 1)
        P[t] = X[t] @ b
    with np.errstate(invalid="ignore"):   # switched off: P = 0, IC 0/0
        v = ic_diario(P[t0:], Y[t0:])
    return pd.Series(np.nan_to_num(v), index=fechas[t0:])

for espera in [21, 63, 126]:
    v = tras_alarma(espera)
    print(f"switched off {espera:3d} sessions: IC 2017-2025 "
          f"{v.mean():.3f} · 2019-S2 to 2020 "
          f"{v['2019-07':'2020-12'].mean():.3f}")

CUSUM alarm (h = 2): 2019-07-19
switched off  21 sessions: IC 2017-2025 0.103 · 2019-S2 to 2020 0.075
switched off  63 sessions: IC 2017-2025 0.101 · 2019-S2 to 2020 0.063
switched off 126 sessions: IC 2017-2025 0.100 · 2019-S2 to 2020 0.054

With 21 sessions switched off, the IC over the 18 months after the change is 0.075, better than any fixed window and 0.016 from the oracle. Waiting longer is costly: with 126 sessions switched off it falls to 0.054, because switching off also switches off f2 and f3, which were still working. Switching off is an emergency measure (not trading with a model under suspicion), not a solution: the value lies in coming back soon with a model estimated on data from the new regime.

movil = {fechas[t].date(): ols(t - 252, t) * 1e4
         for t in sorted(reentreno)}
movil = pd.DataFrame(movil, index=F).T
verdad = b_verdad.set_axis(F, axis=1).set_index(fechas) * 100
fechas_ver = ["2019-07-01", "2019-10-01", "2020-01-02",
              "2022-07-01", "2023-01-03", "2023-07-03", "2024-07-01"]
comp = pd.DataFrame({
    f: {"b1 rolling": movil.loc[pd.Timestamp(f).date(), "f1"],
        "b1 true": verdad.loc[f, "f1"],
        "b2 rolling": movil.loc[pd.Timestamp(f).date(), "f2"],
        "b2 true": verdad.loc[f, "f2"],
        "b5 rolling": movil.loc[pd.Timestamp(f).date(), "f5"],
        "b5 true": verdad.loc[f, "f5"]} for f in fechas_ver}).T
comp.round(1)

            b1 rolling  b1 true  b2 rolling  b2 true  b5 rolling  b5 true
2019-07-01         7.9     -8.0         4.0      5.0         0.9      0.0
2019-10-01         2.4     -8.0         4.9      5.0         1.5      0.0
2020-01-02        -0.5     -8.0         4.7      5.0         0.9      0.0
2022-07-01        -7.2     -8.0         2.6      5.0         1.3      0.0
2023-01-03        -7.9     -8.0         2.2      2.4         3.0      3.1
2023-07-03        -8.2     -8.0         1.8      0.0         4.3      6.0
2024-07-01       -10.1     -8.0        -1.0      0.0         6.6      6.0

The table shows the limits of the 252-session rolling window against the true coefficients. In July 2019 the model still believes in the old effect (7.9 against −8.0), in October it is at 2.4 and in January 2020, at −0.5. In the gradual drift, the coefficient of f2 falls (2.6 in July 2022, 1.8 in July 2023, when the truth is already 0) and that of f5 rises (from 1.3 to 4.3, against 0 and 6): the one-year window runs about six months behind, which is what averaging a year of noisy data costs.

Exercises

1. The 2017 data drift is a textbook covariate shift. The classic correction is to reweight the training data by the density ratio $P_{\text{current}}(x)/P_{\text{base}}(x)$, which is estimated with the adversarial classifier as $p/(1-p)$. Estimate the 2014-2016 OLS with those weights (current = 2017-2019-S1) and compare coefficients and IC with the unweighted fixed model.

actual = (fechas >= "2017-01-01") & (fechas < "2019-07-01")
A = np.vstack([X[base].reshape(-1, 6), X[actual].reshape(-1, 6)])
c = np.r_[np.zeros(base.sum() * NA), np.ones(actual.sum() * NA)]
clf = lgb.LGBMClassifier(
    n_estimators=60, num_leaves=8, min_child_samples=100,
    random_state=0, n_jobs=4, deterministic=True,
    force_row_wise=True, verbose=-1).fit(A, c)
pb = clf.predict_proba(X[base].reshape(-1, 6))[:, 1]
w = pb / (1 - pb)
w = w / w.mean()
sw = np.sqrt(w)
beta_w = np.linalg.lstsq(X1 * sw[:, None], Y[base].ravel() * sw,
                         rcond=None)[0]
print("unweighted (bp):", np.round(beta[1:] * 1e4, 1))
print("weighted (bp):", np.round(beta_w[1:] * 1e4, 1))
print(f"effective size: {w.sum()**2 / (w**2).sum() / w.size:.0%} "
      "of the rows")
for nombre, b in [("unweighted", beta), ("weighted", beta_w)]:
    v = ic_diario(X[actual] @ b[1:], Y[actual]).mean()
    print(f"IC 2017-2019-S1 {nombre}: {v:.3f}")

unweighted (bp): [ 8.9  6.1 -4.5 -0.  -0.   1. ]
weighted (bp): [ 5.5  4.4 -4.3  0.9 -2.   1.2]
effective size: 13% of the rows
IC 2017-2019-S1 unweighted: 0.117
IC 2017-2019-S1 weighted: 0.112

2. The CUSUM monitored the fixed model. Now monitor the 252-session rolling window model with the same detector ($k = 0.06$, $h = 2$, $\mu_0$ = its mean IC in 2017-2018). When does it alarm, and what happens in the gradual drift?

ic_movil = pd.Series(ic_diario(pred["rolling 252"][t0:], Y[t0:]),
                     index=fechas[t0:])
mu_m = ic_movil[:"2018-12-31"].mean()
al = cusum(ic_movil.to_numpy(), mu_m, 0.06, 2.0)
print(f"mean IC 2017-2018: {mu_m:.3f}")
print("alarms:", [str(f.date()) for f in ic_movil.index[al]])
print(ic_movil.groupby(ic_movil.index.year).mean().round(3)
      .to_string())

mean IC 2017-2018: 0.112
alarms: ['2019-07-19', '2019-08-23', '2019-11-26']
2017    0.112
2018    0.112
2019    0.058
2020    0.084
2021    0.097
2022    0.105
2023    0.097
2024    0.133
2025    0.115

3. Design the monitoring of a weekly stock selection model whose label is the 21-session return (overlapping, like the panel's ret_sig_21d). What do you monitor every day, what every week, and how do you calibrate the thresholds?