Unsupervised learning: PCA, clustering and autoencoders.

Structure without a label

The previous lessons in the module learn a function from the features to the label. Unsupervised learning has no label: it looks for structure in the data themselves. In trading it has four common uses, and none of them predicts returns on its own:

flowchart TB
  A[N returns] --> B[PCA]
  A --> C[Groups of assets]
  A --> D[Groups of dates]
  A --> E[Autoencoder]
  B --> F[Factors, residuals]
  C --> G[Diversification]
  D --> H[Regimes]
  E --> I[Anomalies]

• Statistical factors (PCA): summarize 400 series in a few common directions; the residuals are the basis of statistical arbitrage (Avellaneda and Lee, 2010) and the covariance matrix is cleaned with them (lesson 4.1.2).
• Groups of assets: “real” sectors, not those of the classification, to diversify or neutralize a portfolio (HRP, lesson 4.1.4).
• Regimes: group dates with similar behavior (lesson 2.1.7 with hidden Markov models).
• Compression: represent each day with a few variables and measure what does not fit (anomalies).

The problem with evaluating unsupervised learning is that there is no metric like the IC. Here there is, because the data are synthetic and the model that generated them is published (README of public/curso/datos/nivel-5/): the daily return of each of the 400 stocks is

a market factor with one beta per stock, one factor per sector (10 sectors of 40 stocks), idiosyncratic noise with fat tails and small effects (weekly alpha, earnings). Each method is judged by what it recovers of that truth.

The data: 400 daily returns

We work with the simple daily returns of the closes, standardized with the mean and standard deviation of 2011-2018, the fitting period. 2019-2025 is kept to measure out of sample. Standardizing with the whole sample would be the leakage of lesson 5.2.4 in its unsupervised version.

import numpy as np
import pandas as pd

pr = pd.read_parquet("panel5_precios.parquet",
                     columns=["fecha", "ticker", "cierre"])
R = (pr.pivot(index="fecha", columns="ticker", values="cierre")
     .pct_change().iloc[1:])
activos = pd.read_csv("panel5_activos.csv", index_col="ticker")
verdad = pd.read_csv("panel5_activos_verdad.csv", index_col="ticker")
sector = activos.loc[R.columns, "sector"].to_numpy()
beta = verdad.loc[R.columns, "beta"].to_numpy()

ent = (R.index < "2019-01-01")             # 2011-2018
mu, sd = R[ent].mean(), R[ent].std()       # training only
Z = ((R - mu) / sd).to_numpy()
Ze, Zp = Z[ent], Z[~ent]
print("sessions x stocks:", R.shape)
print("training:", Ze.shape[0], " test:", Zp.shape[0])
print("sectors:", pd.Series(sector).value_counts().unique())

sessions x stocks: (3771, 400)
training: 2011  test: 1760
sectors: [40]

PCA and the Marchenko-Pastur edge

PCA looks for the directions of maximum variance: the eigenvectors of the correlation matrix, sorted by their eigenvalue. The practical question is how many components are structure and how many are noise. If returns were independent noise, the eigenvalues of their correlation matrix with $T$ sessions and $N$ assets would fall, for large $T$ and $N$, below

the upper edge of the Marchenko-Pastur distribution. Lesson 4.1.2 used it to clean covariances; here it serves to count factors.

C = np.corrcoef(Ze, rowvar=False)
val, vec = np.linalg.eigh(C)
val, vec = val[::-1], vec[:, ::-1]           # largest to smallest
T, N = Ze.shape
borde = (1 + np.sqrt(N / T)) ** 2            # λ+ with variance 1
print(f"N/T = {N / T:.3f}, Marchenko-Pastur edge = {borde:.2f}")
print("first 12 eigenvalues:", val[:12].round(2))
print(f"above the edge: {(val > borde).sum()}")
print(f"explained variance: PC1 {val[0] / N:.1%}, "
      f"PC1-PC10 {val[:10].sum() / N:.1%}")

N/T = 0.199, Marchenko-Pastur edge = 2.09
first 12 eigenvalues: [120.51   9.79   7.4    6.62   5.88   5.74   5.44   4.71   4.28   3.72
   1.37   1.35]
above the edge: 10
explained variance: PC1 30.1%, PC1-PC10 43.5%

The header figure plots the first twenty. The first, 120.5, takes 30.1% of the variance: it is the market. From the second to the tenth they fall from 9.8 to 3.7 and stay clearly above the edge (2.09). From the eleventh onward they pile up between 1.2 and 1.4: noise. Ten components, which together explain 43.5% of the variance.

The generator has eleven factors. The missing one has not been lost: the market moves all stocks in the same direction, and so does the average of the ten sector factors, so PCA merges those two directions into the first component and what remains distinguishable are nine contrasts between sectors. A principal component is a direction of variance, not a factor of the model one by one; naming each one requires looking at its loadings.

What each component is

pc1 = vec[:, 0] * np.sign(vec[:, 0].sum())
carga = pc1 * sd.to_numpy()          # in return units
f1 = Ze @ pc1                        # component series
ew = R[ent].mean(axis=1).to_numpy()  # equal-weighted portfolio
print(f"corr(PC1 loading, true beta)           = "
      f"{np.corrcoef(carga, beta)[0, 1]:.3f}")
print(f"corr(PC1 loading not rescaled, beta)   = "
      f"{np.corrcoef(pc1, beta)[0, 1]:.3f}")
print(f"corr(PC1 series, equal-weighted port.) = "
      f"{np.corrcoef(f1, ew)[0, 1]:.3f}")

corr(PC1 loading, true beta)           = 0.987
corr(PC1 loading not rescaled, beta)   = 0.620
corr(PC1 series, equal-weighted port.) = 0.999

The PC1 series has a correlation of 0.999 with the equal-weighted portfolio. Its loadings, once multiplied by each stock's volatility (PCA of correlations works with standardized returns), have a correlation of 0.987 with the true beta; without rescaling, only 0.620, because the loading of a correlation measures $\beta_i \sigma_m/\sigma_i$ and the most volatile stocks are penalized. It is a detail often forgotten when using PC1 as a statistical beta.

cargas = pd.DataFrame(vec[:, 1:10] * np.sqrt(val[1:10]),
                      columns=[f"PC{i}" for i in range(2, 11)])
cargas["sector"] = sector
media = cargas.groupby("sector").mean()
print(media.round(2).to_string())

                     PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10
sector
consumo_basico     -0.03 -0.04 -0.04 -0.06 -0.01  0.03 -0.25 -0.16 -0.02
consumo_ciclico    -0.05 -0.03 -0.03  0.03 -0.01  0.01 -0.06  0.19 -0.20
energia             0.46  0.03  0.00 -0.00  0.02 -0.01  0.01 -0.01 -0.01
finanzas           -0.10  0.35  0.10 -0.06  0.05 -0.00  0.04 -0.03 -0.02
industria          -0.04  0.00 -0.03 -0.00  0.01 -0.04 -0.09  0.16  0.21
materiales         -0.06 -0.02 -0.05  0.30  0.03 -0.13  0.05 -0.08  0.00
salud              -0.03 -0.18  0.34 -0.02  0.02  0.00  0.03 -0.02  0.01
servicios_publicos -0.04 -0.05 -0.06 -0.15 -0.25 -0.16  0.09 -0.03 -0.01
tecnologia         -0.03 -0.03 -0.06  0.05 -0.12  0.29  0.10 -0.03  0.03
telecomunicaciones -0.06 -0.13 -0.14 -0.14  0.24 -0.01  0.11 -0.03 -0.01

Each of components 2 to 10 loads heavily on one or two sectors and little on the rest: PC2 is energy (energia, 0.46) against all the others; PC3, financials (finanzas, 0.35) against health care and telecommunications; PC4, health care (salud, 0.34). The last ones mix more (PC8 opposes consumer staples, −0.25, to telecommunications and technology). PCA knows nothing about sectors; it finds them because they are the structure of the covariance. In real data, loadings mix sectors, styles (value, momentum) and size, and their interpretation changes from one window to another (lesson 4.1.6).

Grouping stocks

If each sector occupies its own direction in the space of components 2 to 10, a k-means on those loadings (normalized so that the direction counts, not the size) should find the groups. To measure it, the adjusted Rand index (ARI): 1 if two partitions coincide (with any group names), around 0 if they are as alike as they would be by chance.

from sklearn.cluster import KMeans
from sklearn.metrics import adjusted_rand_score

L = vec[:, 1:10] * np.sqrt(val[1:10])
L = L / np.linalg.norm(L, axis=1, keepdims=True)
km = KMeans(n_clusters=10, n_init=10, random_state=0).fit(L)
print(f"ARI k-means (PC2-PC10 loadings) versus sectors: "
      f"{adjusted_rand_score(sector, km.labels_):.3f}")
print(pd.crosstab(km.labels_, sector).max(axis=1).tolist())

ARI k-means (PC2-PC10 loadings) versus sectors: 1.000
[40, 40, 40, 40, 40, 40, 40, 40, 40, 40]

ARI of 1.000: the ten groups are the ten sectors, with 40 stocks each. The most common alternative in finance is hierarchical clustering on the correlation distance $d_{ij} = \sqrt{2(1-\rho_{ij})}$, the same one used by the HRP of lesson 4.1.4. The cell tries it with and without the first component:

from sklearn.cluster import AgglomerativeClustering

def jerarquico(Zv, quitar_mercado=True, k=10):
    """Clusters by correlation (distance sqrt(2(1-rho)))."""
    if quitar_mercado:
        w, v = np.linalg.eigh(np.corrcoef(Zv, rowvar=False))
        Zv = Zv - np.outer(Zv @ v[:, -1], v[:, -1])
    rho = np.corrcoef(Zv, rowvar=False)
    dist = np.sqrt(np.clip(2 * (1 - rho), 0, None))
    return AgglomerativeClustering(
        n_clusters=k, metric="precomputed",
        linkage="average").fit(dist).labels_

for quitar in (False, True):
    lab = jerarquico(Ze, quitar)
    print(f"remove PC1={quitar!s:5}: ARI "
          f"{adjusted_rand_score(sector, lab):.3f}, sizes "
          f"{sorted(np.bincount(lab).tolist(), reverse=True)}")

remove PC1=False: ARI 0.001, sizes [386, 5, 2, 1, 1, 1, 1, 1, 1, 1]
remove PC1=True : ARI 1.000, sizes [40, 40, 40, 40, 40, 40, 40, 40, 40, 40]

Without removing anything, the result is useless: one group of 386 stocks and nine of one to five (ARI 0.001). With an average correlation of 0.32 because of the market, all stocks look alike and the sector differences stay buried; average linkage keeps merging everything into a single group. Removing the projection on PC1, the residuals keep only the sector part and the clustering recovers the ten exact sectors.

That it comes out perfect says more about the data than about the method: the synthetic sectors are clean (40 stocks, one factor each, no company straddling two) and there are eight years of daily data. The next cell repeats the clustering with shorter windows, in six different stretches of 2011-2018:

filas = {}
for n in (21, 63, 126, 252, 504):
    aris = [adjusted_rand_score(sector, jerarquico(Ze[s:s + n]))
            for s in range(0, 1800, 300)]
    filas[n] = {"ARI medio": np.mean(aris), "mínimo": np.min(aris),
                "máximo": np.max(aris)}
pd.DataFrame(filas).T.rename_axis("sesiones").round(2)

          ARI medio  mínimo  máximo
sesiones
21             0.21    0.11    0.29
63             0.39    0.29    0.55
126            0.58    0.46    0.82
252            0.84    0.74    0.92
504            0.97    0.93    0.99

With one month of data (21 sessions), the average ARI is 0.21; with a quarter, 0.39; with a year, 0.84; with two, 0.97. The groups of a clustering done with little data depend on the window, and an HRP recomputed with 63 sessions reorganizes its tree every time. In real data, with a less clean structure that also changes, the window length is a trade-off between noise and recency.

The first component in crises

The variance taken by PC1 is not constant. Computed quarter by quarter, it measures how much stocks move together: it is the idea of the absorption ratio of Kritzman et al. (2011) as a measure of systemic risk.

Rv = R.to_numpy()
trim = {}
for per, idx in R.groupby(R.index.to_period("Q")).indices.items():
    Zq = Rv[idx]
    Zq = (Zq - Zq.mean(0)) / Zq.std(0)
    ev = np.linalg.eigvalsh(np.corrcoef(Zq, rowvar=False))
    trim[str(per)] = ev[-1] / ev.sum()
trim = pd.Series(trim)
print(f"PC1 by quarter: median {trim.median():.2f}, "
      f"minimum {trim.min():.2f} ({trim.idxmin()}), "
      f"maximum {trim.max():.2f} ({trim.idxmax()})")
print(trim.nlargest(5).round(2).to_string())

PC1 by quarter: median 0.31, minimum 0.13 (2019Q2), maximum 0.68 (2020Q1)
2020Q1    0.68
2020Q2    0.68
2015Q4    0.53
2018Q4    0.53
2016Q1    0.52

The median is 31%, with a minimum of 13% in the second quarter of 2019 and a maximum of 68% in the first two quarters of 2020, the generator's hardest stress stretch (market volatility of 80% annualized). The other peaks (late 2015, late 2018, 2022) coincide with the other stress stretches. It is what you expect of any portfolio in a crisis: diversification across stocks disappears just when it is needed (lesson 4.2.4, correlation stress testing). Note that quarters of 63 sessions give a Marchenko-Pastur edge of about 13: with short windows only the market stands out from the noise (exercise 2 counts it year by year).

Regimes with clustering

What groups stocks can also group dates. Each Friday is described with four features computed only with the past: the logarithm of market volatility over 21 sessions, the approximate average correlation (market variance divided by the average variance of the stocks, 63 sessions), the market return over 21 sessions and the logarithm of the cross-sectional dispersion between stocks. A mixture of three Gaussians is fitted with 2012-2018 and classifies all the weeks; the generator publishes the true regime of each one (calma for calm, estres for stress, rebote for rebound) in panel5_coef_verdad.csv.

from sklearn.mixture import GaussianMixture
from sklearn.preprocessing import StandardScaler

reg = pd.read_csv("panel5_coef_verdad.csv", parse_dates=["fecha"],
                  index_col="fecha")["regimen"]
mkt = R.mean(axis=1)
var_i = R.rolling(63).var().mean(axis=1)
F = pd.DataFrame({
    "log_vol_21": np.log(mkt.rolling(21).std() * np.sqrt(252)),
    "correlacion_63": mkt.rolling(63).var() / var_i,
    "ret_21": mkt.rolling(21).sum(),
    "log_dispersion_21": np.log(R.std(axis=1).rolling(21).mean()),
}).reindex(reg.index, method="ffill")    # Friday of each week
e2 = F.index < "2019-01-01"
esc = StandardScaler().fit(F[e2])
gm = GaussianMixture(3, covariance_type="full", n_init=5,
                     random_state=0).fit(esc.transform(F[e2]))
grupo = pd.Series(gm.predict(esc.transform(F)), index=F.index)
print(pd.DataFrame(esc.inverse_transform(gm.means_),
                   columns=F.columns).round(3))
print(pd.crosstab(reg[~e2], grupo[~e2],
                  rownames=["2019-2025"], colnames=["grupo"]))

   log_vol_21  correlacion_63  ret_21  log_dispersion_21
0      -1.691           0.324   0.035             -4.070
1      -2.224           0.186   0.008             -4.158
2      -1.242           0.425  -0.052             -3.936
grupo       0    1   2
2019-2025
calma      72  190   1
estres     15    4  33
rebote     36    0  15

The group numbers are assigned by the algorithm; their means say what each one is. Group 2 has the highest volatility ($e^{-1.24} \approx$ 29% annualized), the highest correlation (0.43) and a negative return: it is stress. Group 1 is quiet calm (volatility of 11%, correlation 0.19). Group 0 sits between the two. In 2019-2025, outside the fitting sample:

• Stress: 33 of its 52 weeks fall in group 2 and 15 in the intermediate one.
• Calm: 190 of 263 in group 1 and 72 in the intermediate one; a single one in the stress group.
• Rebound: none in calm, but spread between the intermediate group (36) and stress (15). A rebound is volatile like stress and is only distinguished by the sign of the return, which over three months of data is noisy.

Grouping dates also has its version of delay: the features are backward-looking windows, so the regime has to last to move them.

estres = grupo == int(np.argmax(gm.means_[:, 0]))
inicios = reg.index[(reg == "estres") & (reg.shift() != "estres")]
for ini in inicios[inicios >= "2019-01-01"]:
    tras = estres[ini:]
    primera = tras.index[tras.to_numpy()][0]
    print(f"stress from {ini:%Y-%m-%d}: the group flags it on "
          f"{primera:%Y-%m-%d} ({(primera - ini).days // 7} weeks)")

stress from 2020-02-21: the group flags it on 2020-03-06 (2 weeks)
stress from 2022-01-07: the group flags it on 2022-02-11 (5 weeks)
stress from 2025-02-21: the group flags it on 2025-03-14 (3 weeks)

The February 2020 stress appears in group 2 two weeks after it starts; that of January 2022, five (it was a mild stress, 25% volatility); that of February 2025, three. A detector without future information always arrives late: the useful question is how late and how much that delay costs the strategy that uses it. Lesson 2.1.7 did the same with an HMM, which adds the persistence of the states, and lesson 5.10.2 returns to regimes as data drift.

Autoencoders

An autoencoder is a network that learns to reproduce its input by passing through a narrow bottleneck: an encoder compresses the day's 400 returns into $k$ numbers and a decoder reconstructs them. It is trained with squared error, without a label.

flowchart TB
  A[400 returns] --> B[tanh layer 64]
  B --> C[Bottleneck of 10]
  C --> D[tanh layer 64]
  D --> E[400 reconstructed]
  E --> F[Squared error]

If all the layers are linear, the best autoencoder with a bottleneck of $k$ reconstructs from the same subspace as the first $k$ principal components: it is PCA with another algorithm. The promise of nonlinear layers is to capture structure that PCA does not see. The cell trains both in Keras 3 with the PyTorch backend (the KERAS_BACKEND variable is set by the runner; lesson 5.6.1 explains the backends), with a seed (keras.utils.set_random_seed), deterministic torch algorithms and 60 epochs on 2011-2018. In production you would try more architectures, regularization and early stopping with a purged validation.

import os
import keras
import torch

torch.use_deterministic_algorithms(True)
print("backend:", keras.backend.backend(), os.environ["KERAS_BACKEND"])

def r2(A, B):
    """Fraction of the (standardized) variance reconstructed."""
    return 1 - ((A - B) ** 2).sum() / (A ** 2).sum()

def autoencoder(k, oculta=None, epocas=60):
    keras.utils.set_random_seed(0)
    x = entrada = keras.Input((N,))
    if oculta:
        x = keras.layers.Dense(oculta, activation="tanh")(x)
    x = keras.layers.Dense(k, name="cuello")(x)
    if oculta:
        x = keras.layers.Dense(oculta, activation="tanh")(x)
    salida = keras.layers.Dense(N)(x)
    m = keras.Model(entrada, salida)
    m.compile(optimizer=keras.optimizers.Adam(1e-3), loss="mse")
    m.fit(Ze.astype("float32"), Ze.astype("float32"), epochs=epocas,
          batch_size=64, shuffle=True, verbose=0)
    return m

k = 10
V = vec[:, :k]
res = {"PCA": (r2(Ze, Ze @ V @ V.T), r2(Zp, Zp @ V @ V.T))}
for nombre, oculta in [("autoencoder lineal", None),
                       ("autoencoder tanh 64", 64)]:
    m = autoencoder(k, oculta)
    pe = m.predict(Ze.astype("float32"), verbose=0)
    pp = m.predict(Zp.astype("float32"), verbose=0)
    res[nombre] = (r2(Ze, pe), r2(Zp, pp))
pd.DataFrame(res, index=["R2 2011-2018", "R2 2019-2025"]).T.round(4)

backend: torch torch
                     R2 2011-2018  R2 2019-2025
PCA                        0.4352        0.4811
autoencoder lineal         0.4350        0.4811
autoencoder tanh 64        0.4301        0.4408

The measure is the fraction of the standardized variance that is reconstructed ($R^2$), in sample (2011-2018) and out of sample (2019-2025):

• The linear autoencoder matches PCA to the fourth decimal out of sample (0.4811 both): it learns the same subspace, more slowly.
• The nonlinear one (tanh layers of 64) loses little in sample (0.4301 versus 0.4352) and quite a lot out of sample (0.4408 versus 0.4811). It has nothing nonlinear to find, because the true model is a linear factor model, and its extra parameters fit peculiarities of 2011-2018 that do not repeat.

curva = {}
for kk in (1, 2, 3, 5, 10, 15, 20, 30, 50):
    Vk = vec[:, :kk]
    curva[kk] = {"dentro": r2(Ze, Ze @ Vk @ Vk.T),
                 "fuera": r2(Zp, Zp @ Vk @ Vk.T)}
pd.DataFrame(curva).T.rename_axis("componentes").round(4)

             dentro   fuera
componentes
1            0.3013  0.3678
2            0.3257  0.3894
3            0.3442  0.4046
5            0.3755  0.4281
10           0.4352  0.4811
15           0.4517  0.4889
20           0.4672  0.4963
30           0.4963  0.5116
50           0.5500  0.5412

The figure places the two autoencoders on the PCA curve. With a bottleneck of 10 (the components above the edge), the nonlinear one reconstructs out of sample little more than PCA with 5 components (0.428). And the PCA curve has its own lesson: going from 10 to 50 components raises the in-sample $R^2$ (from 0.435 to 0.550) more than the out-of-sample one (from 0.481 to 0.541), because those 40 extra components are the noise that the Marchenko-Pastur edge had already discarded.

That the $R^2$ is higher out of sample than in sample is not leakage: 2019-2025 contains 2020, when the market took up to 68% of the variance, and a period with more correlation is reconstructed better with the same components. That is why methods are compared over the same period.

Common mistakes

• PCA on non-standardized returns or with the whole sample. Without standardizing, volatile stocks dominate the components; standardizing with the whole sample, the test data enter the fit.
• Counting factors by eye. The Marchenko-Pastur edge gives a reference; the components below it are noise and adding them raises the in-sample fit, not the out-of-sample one.
• Clustering correlations without removing the market. You get one giant group and a few stragglers.
• Reading group numbers as fixed labels. k-means or a mixture may number the groups differently when refitted; a group is identified by its means.
• Regimes fitted with the whole sample. A detector that “gets” past regimes right with later information does not say when it would have detected them.
• Using an autoencoder because it is more modern. If it does not improve on PCA out of sample with the same dimension, it is PCA that is more expensive and less stable.

Exercises

1. Repeat the hierarchical clustering without removing the market but with 11 groups (in case the market takes one). Does it improve? Then try Ward linkage on the residuals without PC1.

lab = jerarquico(Ze, quitar_mercado=False, k=11)
print(sorted(np.bincount(lab).tolist(), reverse=True))
print(f"ARI: {adjusted_rand_score(sector, lab):.3f}")
lab_w = AgglomerativeClustering(n_clusters=10, linkage="ward").fit(
    (Ze - np.outer(Ze @ vec[:, 0], vec[:, 0])).T).labels_
print(f"Ward on the residuals: ARI "
      f"{adjusted_rand_score(sector, lab_w):.3f}")

[385, 5, 2, 1, 1, 1, 1, 1, 1, 1, 1]
ARI: 0.001
Ward on the residuals: ARI 0.977

2. Compute year by year how many eigenvalues exceed the Marchenko-Pastur edge and what fraction of the variance PC1 takes. What happens in 2020?

anios = {}
for a, idx in R.groupby(R.index.year).indices.items():
    Za = Rv[idx]
    Za = (Za - Za.mean(0)) / Za.std(0)
    ev = np.linalg.eigvalsh(np.corrcoef(Za, rowvar=False))
    b = (1 + np.sqrt(N / len(idx))) ** 2
    anios[a] = {"sesiones": len(idx), "borde": round(b, 2),
                "encima": int((ev > b).sum()),
                "PC1": round(ev[-1] / N, 2)}
pd.DataFrame(anios).T.astype({"sesiones": int, "encima": int})

      sesiones  borde  encima   PC1
2011       251   5.12       8  0.31
2012       250   5.13       7  0.32
2013       252   5.11       8  0.23
2014       252   5.11       7  0.28
2015       252   5.11       8  0.38
2016       252   5.11       7  0.37
2017       251   5.12       8  0.20
2018       251   5.12       8  0.34
2019       252   5.11       8  0.26
2020       253   5.10       2  0.58
2021       252   5.11       7  0.28
2022       251   5.12       6  0.44
2023       250   5.13       9  0.23
2024       252   5.11       8  0.29
2025       250   5.13       7  0.35

3. A colleague computes the first five principal components with the returns of 2011-2025, uses each stock's exposure to each component as a feature of a 5.5.3 model and validates it with walk-forward from 2019. What is wrong?