Time series architectures: TCN, N-BEATS, PatchTST and iTransformer.

The test bench: volatility with a known answer

Lesson 5.6.3 trained 1D CNNs and LSTMs on windows. Since 2018, architectures designed for time series have appeared that are cited in almost every forecasting paper: the TCN (dilated causal convolutions), N-BEATS and its successor N-HiTS (dense blocks with double residuals) and transformers adapted to time series, such as PatchTST and iTransformer. This lesson codes small versions of the four, trains N-HiTS with a library and compares them where they can really be measured.

The problem is the one from lesson 5.7.1: forecast the logarithm of tomorrow's realized variance (log RV) for the six assets of vol5_diario.parquet. Now each model receives the window of the last 64 days of each series, not the precomputed HAR features, and has to find the structure by itself. The data are synthetic and the generator is a HAR on the logarithm of integrated variance, so there are three benchmarks:

import copy
import logging
import warnings
import numpy as np
import pandas as pd
import torch
from torch import nn

torch.manual_seed(0)
torch.use_deterministic_algorithms(True)

v = pd.read_parquet("vol5_diario.parquet")
P = v.pivot(index="fecha", columns="activo", values="rv_5min")
P = np.log(P)
ACT = list(P.columns)
fechas = P.index
en_ent = fechas < "2016-01-01"
mu, sd = P[en_ent].mean(), P[en_ent].std()   # training only
Z = ((P - mu) / sd).to_numpy(dtype=np.float32)   # (days, 6)

L = 64
pos = np.arange(L, len(fechas))     # position of the target day
f_obj = fechas[pos]
tramo = np.select([f_obj < "2016-01-01", f_obj < "2018-01-01"],
                  ["ent", "val"], "prueba")

def ventanas(ts):
    W = np.stack([Z[t - L:t] for t in ts])     # (n, L, 6)
    return torch.tensor(W), torch.tensor(Z[ts])

# multivariate: one sample per date (for iTransformer)
D = {k: ventanas(pos[tramo == k]) for k in ("ent", "val", "prueba")}
# univariate: one sample per date and asset (the rest)
U = {k: (W.permute(0, 2, 1).reshape(-1, L), Y.reshape(-1))
     for k, (W, Y) in D.items()}
print({k: tuple(W.shape) for k, (W, _) in U.items()})

def har_x(W):
    return np.c_[np.ones(len(W)), W[:, -1], W[:, -5:].mean(1),
                 W[:, -22:].mean(1)].astype(float)

b_har = np.linalg.lstsq(har_x(U["ent"][0].numpy()),
                        U["ent"][1].double().numpy(), rcond=None)[0]

ver = pd.read_parquet("vol5_verdad.parquet")
log_iv = np.log(ver.pivot(index="fecha", columns="activo",
                          values="iv_prevista"))[ACT]
cte = (P - log_iv)[en_ent].mean()   # RV versus IV: a constant
ORA = ((log_iv + cte - mu) / sd).to_numpy()

y_p = U["prueba"][1].double().numpy()
ref = {"naive": U["prueba"][0][:, -1].double().numpy(),
       "HAR": har_x(U["prueba"][0].numpy()) @ b_har,
       "oracle": ORA[pos[tramo == "prueba"]].reshape(-1)}
ecm_ref = {k: np.mean((y_p - p) ** 2) for k, p in ref.items()}
print("HAR:", b_har.round(3))
print({k: round(e, 3) for k, e in ecm_ref.items()})

{'ent': (11706, 64), 'val': (3018, 64), 'prueba': (12066, 64)}
HAR: [0.001 0.241 0.341 0.308]
{'naive': np.float64(0.842), 'HAR': np.float64(0.577), 'oracle': np.float64(0.532)}

• Each series is standardized with the mean and standard deviation of its own training period (2008-2015): that way the six share a scale and a global model can learn from all of them.
• The splits go by target date, as in 5.7.1 (ent, val and prueba are training, validation and test). With a one-day horizon there is no need to purge: no training window has its target in the validation period.
• The two views of the same data matter: 11,706 univariate windows (date × asset) for the models that process one series at a time, and 1,951 multivariate windows (one per date, six series) for the iTransformer.
• The oracle needs a constant per asset (RV does not see the overnight part and its logarithm is biased), estimated on training. With it, its test error is the floor of the table.

The test error (in squared standardized units) is 0.842 for the naive forecast, 0.577 for the HAR and 0.532 for the oracle. The gap a network can close is the one between the HAR and the oracle, 0.045: 8% of the HAR's error. Everything else is RV measurement noise and shocks nobody can forecast.

TCN: dilated causal convolutions

A temporal convolutional network (TCN; Bai, Kolter and Koltun, 2018) stacks 1D convolutions with two conditions:

• Causal: the output on day $t$ only uses inputs from $t$ backward. It is achieved by zero-padding on the left only.
• Dilated: layer $j$ looks at positions $d_j$ days apart ($d = 1, 2, 4, \ldots$). The receptive field grows exponentially with the number of layers:

The figure draws it with three layers (dilations 1, 2 and 4): 8 days. With kernel $k = 2$ and dilations 1, 2, 4, 8, 16 and 32, $R = 64$: six layers are enough to see the whole window. Each layer adds its output to its input (residual connection), which makes deep networks easier to train. The network has no recurrence: it processes all positions in parallel, and that is why it trains faster than an LSTM.

class TCN(nn.Module):
    def __init__(self, canales=16, k=2, dil=(1, 2, 4, 8, 16, 32)):
        super().__init__()
        self.k, self.dil = k, dil
        self.entrada = nn.Conv1d(1, canales, 1)
        self.capas = nn.ModuleList(
            [nn.Conv1d(canales, canales, k, dilation=d) for d in dil])
        self.salida = nn.Linear(canales, 1)

    def forward(self, x):                      # x: (n, L)
        h = self.entrada(x.unsqueeze(1))       # (n, channels, L)
        for conv, d in zip(self.capas, self.dil):
            relleno = nn.functional.pad(h, ((self.k - 1) * d, 0))
            h = h + torch.relu(conv(relleno))  # causal and residual
        return self.salida(h[:, :, -1]).squeeze(-1)

torch.manual_seed(0)
tcn = TCN()
x = torch.zeros(1, 96, requires_grad=True)     # longer than 64
tcn(x).sum().backward()
vistas = torch.nonzero(x.grad[0]).flatten()
print("positions with gradient:", len(vistas),
      f"(from {vistas.min().item()} to {vistas.max().item()})")
print("parameters:", sum(p.numel() for p in tcn.parameters()))

positions with gradient: 64 (from 32 to 95)
parameters: 3217

With a 96-day input, the output only has a gradient with respect to the last 64 positions: the receptive field is exactly the one in the formula, and nothing from the future gets in (there are no positions after the last one). It is the autograd check of lesson 5.7.1 used to audit an architecture.

N-BEATS and N-HiTS: blocks that subtract what they explain

N-BEATS (Oreshkin et al., 2020) uses neither convolutions nor recurrence: only dense layers organized in blocks. Each block receives a window and produces two things: a backcast (its reconstruction of the window) and a forecast. The next block receives the window minus the backcast, that is, what the previous one has not explained, and the final forecast is the sum of the forecasts of all the blocks:

flowchart TB
  X[Window x] --> B1[Block 1]
  B1 -- backcast --> R1[x minus backcast]
  B1 -- forecast --> S[Sum]
  R1 --> B2[Block 2]
  B2 -- forecast --> S
  B2 -- backcast --> R2[residual]
  R2 --> B3[Block 3]
  B3 -- forecast --> S

In the “interpretable” version, the blocks are restricted to trend bases (polynomials) and seasonality bases (sines and cosines), and their forecasts can be read separately. N-HiTS (Challu et al., 2023) adds two ideas for long horizons: each block subsamples the input at a different frequency (pooling) and forecasts a few points that are then interpolated, so some blocks specialize in the slow parts and others in the fast ones, with much less computation. With a one-day horizon those advantages do not appear; the neuralforecast section checks it.

class Bloque(nn.Module):
    def __init__(self, L, H, ancho):
        super().__init__()
        self.mlp = nn.Sequential(nn.Linear(L, ancho), nn.ReLU(),
                                 nn.Linear(ancho, ancho), nn.ReLU())
        self.backcast = nn.Linear(ancho, L)
        self.forecast = nn.Linear(ancho, H)

    def forward(self, x):
        h = self.mlp(x)
        return self.backcast(h), self.forecast(h)

class NBeats(nn.Module):
    def __init__(self, L=64, H=1, ancho=64, n_bloques=3):
        super().__init__()
        self.bloques = nn.ModuleList(
            [Bloque(L, H, ancho) for _ in range(n_bloques)])

    def forward(self, x):                      # x: (n, L)
        residuo, prevision = x, 0
        for b in self.bloques:
            bc, fc = b(residuo)
            residuo = residuo - bc
            prevision = prevision + fc
        return prevision.squeeze(-1)

torch.manual_seed(0)
print("parameters:", sum(p.numel() for p in NBeats().parameters()))

parameters: 37635

Attention: what each token looks at

A transformer (Vaswani et al., 2017) transforms a sequence of tokens (vectors of dimension $d$) by letting each one look at all the others. Each token produces a query $q$, a key $k$ and a value $v$ (learned linear projections), and a token's output is the average of everyone's values, weighted by how similar its query is to each key:

The $\sqrt{d}$ factor keeps the dot products from growing with the dimension and saturating the softmax. Several heads repeat the computation with different projections, and a positional encoding added to each token tells the network the order, which attention alone does not see. To make it concrete, the tokens below are the eight 8-day patches of the index's last test window, without learned projections ($Q = K = V$):

indice = ACT.index("indice")
fichas = U["prueba"][0].reshape(-1, 6, L)[-1, indice].reshape(8, 8)
Q = K = V = fichas.unsqueeze(0)                # (1, 8 tokens, d=8)
pesos = torch.softmax(Q @ K.transpose(-2, -1) / 8 ** 0.5, dim=-1)
a_mano = pesos @ V
de_torch = nn.functional.scaled_dot_product_attention(Q, K, V)
print("equal:", torch.allclose(a_mano, de_torch, atol=1e-6))
print("weights of the last patch over the eight:")
print(pesos[0, -1].numpy().round(2))

equal: True
weights of the last patch over the eight:
[0.08 0.07 0.14 0.08 0.05 0.08 0.12 0.38]

Each row of weights sums to 1 and says which patches a token takes information from: the last patch puts 0.38 on itself and spreads the rest; after it, the third (0.14) and the seventh (0.12) weigh more, the ones most similar to it. In a real transformer, the learned projections decide what “similar” means. The cost is quadratic in the number of tokens: with one token per day and one-year windows it would be 252² products per layer and head. That cost, and the fact that a single day carries very little information, explain why the first transformers for time series (Informer, Autoformer, FEDformer) performed worse than they promised: Zeng et al. (2023) showed that a linear model on the window (DLinear) beat them on the usual benchmarks. The two answers that follow change what a token is.

PatchTST and iTransformer: what a token is

PatchTST (Nie et al., 2023) cuts each series into patches of several consecutive days and uses each patch as a token: with 64 days and patches of 8, attention works over 8 tokens instead of 64, and each token carries a piece of local dynamics. It is also channel-independent: each series is processed separately with the same weights, like the univariate models of this lesson. The title of its paper, “a time series is worth 64 words”, comes from there.

iTransformer (Liu et al., 2024) turns the idea around: each complete series (the 64 days of one asset) is projected to a single token, and attention relates series to each other. The time dimension is processed by the dense network that the transformer applies to each token. It makes sense when some series carry information about others: the volatility of the index about that of a stock, or that of one exchange rate about another.

flowchart LR
  S[Series 64 days] --> P[8 patches]
  P --> T1[8 tokens]
  T1 --> A1[PatchTST]
  M[6 series] --> E[1 token per series]
  E --> A2[iTransformer]

def codificador(d, cabezas=4, capas=2):
    capa = nn.TransformerEncoderLayer(d, cabezas, dim_feedforward=64,
                                      dropout=0.1, batch_first=True)
    return nn.TransformerEncoder(capa, capas,
                                 enable_nested_tensor=False)

class PatchTST(nn.Module):
    def __init__(self, L=64, parche=8, d=32):
        super().__init__()
        self.parche, n = parche, L // parche
        self.proy = nn.Linear(parche, d)
        self.posicion = nn.Parameter(0.02 * torch.randn(n, d))
        self.enc = codificador(d)
        self.cabeza = nn.Linear(n * d, 1)

    def forward(self, x):                      # x: (n, L)
        fichas = x.unfold(1, self.parche, self.parche)  # (n, 8, 8)
        h = self.enc(self.proy(fichas) + self.posicion)
        return self.cabeza(h.flatten(1)).squeeze(-1)

class ITransformer(nn.Module):
    def __init__(self, L=64, d=32):
        super().__init__()
        self.proy = nn.Linear(L, d)            # one series, one token
        self.enc = codificador(d)
        self.cabeza = nn.Linear(d, 1)

    def forward(self, x):                      # x: (n, L, 6)
        h = self.enc(self.proy(x.transpose(1, 2)))   # (n, 6, d)
        return self.cabeza(h).squeeze(-1)            # (n, 6)

torch.manual_seed(0)
for m in (PatchTST(), ITransformer()):
    print(type(m).__name__, "parameters:",
          sum(p.numel() for p in m.parameters()))

PatchTST parameters: 17889
ITransformer parameters: 19201

Both use PyTorch's standard layers (nn.TransformerEncoderLayer: multi-head attention, dense network, normalization and residuals). enable_nested_tensor=False avoids a warning: the nested-tensor optimization only applies with padding masks, and there are none here.

Training and comparing

The four networks are trained with the same loop as in 5.7.1, with one practical difference: since the data are already tensors in memory, batches are drawn with a dedicated permutation (with its own generator) instead of a DataLoader; the result is the same and it avoids overhead. Early stopping on validation (2016-2017), test on 2018-2025. They are small networks with few epochs so that the lesson runs in a couple of minutes on CPU; in production you would train longer and with several seeds.

def entrenar(modelo, datos, epocas=15, lr=1e-3, paciencia=4,
             lote=256, semilla=0):
    Xe, Ye = datos["ent"]
    Xv, Yv = datos["val"]
    opt = torch.optim.Adam(modelo.parameters(), lr=lr)
    gen = torch.Generator().manual_seed(semilla)
    mejor, sin_mejora = (np.inf, None, -1), 0
    for epoca in range(epocas):
        modelo.train()
        orden = torch.randperm(len(Xe), generator=gen)
        for i in range(0, len(Xe), lote):
            j = orden[i:i + lote]
            opt.zero_grad()
            nn.functional.mse_loss(modelo(Xe[j]), Ye[j]).backward()
            opt.step()
        modelo.eval()
        with torch.no_grad():
            val = nn.functional.mse_loss(modelo(Xv), Yv).item()
        if val < mejor[0]:
            mejor = (val, copy.deepcopy(modelo.state_dict()), epoca)
            sin_mejora = 0
        else:
            sin_mejora += 1
            if sin_mejora >= paciencia:
                break
    modelo.load_state_dict(mejor[1])
    return mejor[2]

def prever(modelo, X):
    modelo.eval()
    with torch.no_grad():
        return modelo(X).reshape(-1).double().numpy()

redes = {"TCN": (TCN, U, {}), "N-BEATS": (NBeats, U, {}),
         "PatchTST": (PatchTST, U, {}),
         "iTransformer": (ITransformer, D,
                          {"epocas": 40, "paciencia": 8, "lote": 32})}
pred, filas = dict(ref), []
for nombre, (clase, datos, kw) in redes.items():
    torch.manual_seed(0)
    m = clase()
    epoca = entrenar(m, datos, **kw)
    pred[nombre] = prever(m, datos["prueba"][0])
    filas.append((nombre, sum(p.numel() for p in m.parameters()),
                  epoca))

hueco = ecm_ref["HAR"] - ecm_ref["oracle"]
tabla = pd.DataFrame(filas, columns=["model", "parameters",
                                     "best epoch"]).set_index("model")
tabla["test MSE"] = [np.mean((y_p - pred[n]) ** 2)
                     for n in tabla.index]
tabla["gap closed"] = (ecm_ref["HAR"] - tabla["test MSE"]) / hueco
tabla.round(3)

              parameters  best epoch  test MSE  gap closed
model
TCN                 3217           7     0.570       0.145
N-BEATS            37635           4     0.590      -0.290
PatchTST           17889           7     0.588      -0.249
iTransformer       19201           4     0.623      -1.043

The iTransformer's predictions come out in date × asset order, the same as the univariate ones (reshape of (n, 6)), so they are all compared row by row against the same targets. The “gap closed” column measures what fraction of the distance between the HAR and the oracle each model covers: 0 is the HAR, 1 would be the oracle and a negative value, worse than the HAR.

semana = np.repeat(f_obj[tramo == "prueba"].to_period("W"), 6)
for nombre in ("TCN", "N-BEATS"):
    dif = (y_p - pred["HAR"]) ** 2 - (y_p - pred[nombre]) ** 2
    s = pd.Series(dif).groupby(semana).mean()
    print(f"{nombre} versus HAR: {dif.mean():+.4f} per day; "
          f"t = {s.mean() / s.std() * np.sqrt(len(s)):.1f}")

TCN versus HAR: +0.0065 per day; t = 3.1
N-BEATS versus HAR: -0.0130 per day; t = -5.7

The same family in neuralforecast

Coding the architecture teaches what it does; for actual work you use a library. neuralforecast (Nixtla) ships more than thirty models (TCN, NBEATS, NHITS, PatchTST, iTransformer, DeepAR…) on top of PyTorch and Lightning, with a single data format: a long table with unique_id (the series), ds (the date or an integer) and y. Its cross_validation with refit=False does exactly the scheme of this lesson: it trains once on everything before validation, stops early on validation and forecasts each test day with the window that ends the previous day.

An integer (the session position) is used as ds instead of the date: NYSE sessions do not have a regular calendar frequency, and freq=1 keeps the library from inventing dates. NHITS's default configuration has about 2.5 million parameters and is designed for long horizons; with 64-neuron layers it drops to about 58,000, trains in a fraction of the time and forecasts better.

from neuralforecast import NeuralForecast
from neuralforecast.models import NHITS

for n in ("pytorch_lightning", "lightning.pytorch",
          "lightning_fabric", "lightning", "neuralforecast"):
    logging.getLogger(n).setLevel(logging.ERROR)
warnings.filterwarnings("ignore", module="pytorch_lightning")
warnings.filterwarnings("ignore", module="lightning")

larga = (pd.DataFrame(Z, columns=ACT).rename_axis("ds").reset_index()
         .melt(id_vars="ds", var_name="unique_id", value_name="y"))
nhits = NHITS(h=1, input_size=L, mlp_units=[[64, 64]] * 3,
              max_steps=600, val_check_steps=50,
              early_stop_patience_steps=4, learning_rate=1e-3,
              batch_size=32, windows_batch_size=256, scaler_type=None,
              random_seed=0, accelerator="cpu", deterministic=True,
              logger=False, enable_progress_bar=False,
              enable_model_summary=False, enable_checkpointing=False)
nf = NeuralForecast(models=[nhits], freq=1)
cv = nf.cross_validation(larga, val_size=(tramo == "val").sum(),
                         test_size=(tramo == "prueba").sum(),
                         n_windows=None, step_size=1, refit=False)
print(cv.shape, "parameters:",
      sum(p.numel() for p in nf.models[0].parameters()))
cv = cv.sort_values(["ds", "unique_id"])       # date × asset
assert np.allclose(cv["y"].to_numpy(), y_p, atol=1e-6)
pred["NHITS (neuralforecast)"] = cv["NHITS"].to_numpy()
print(f"NHITS test MSE: {np.mean((y_p - cv['NHITS']) ** 2):.3f}")

(12066, 5) parameters: 58499
NHITS test MSE: 0.600

The assert checks that the library's targets are the same, in the same order, as those of the previous table: if they were not, the comparison would be worthless. It is a line worth adding whenever results from your own code and from a library are mixed.

What the results say

The figure at the top gathers the test results. Four readings:

• All the networks comfortably beat the naive forecast and none comes close to the oracle. The gap between the HAR and the oracle is small (8% of the error) and it is almost all there is: the rest of the error is measurement and shocks.
• Only the TCN improves on the HAR, and by little: 0.0065 per day, 15% of the gap, with t = 3.1 on the weekly difference. It is moderate evidence from a single training run and one split, after trying five models; N-BEATS and PatchTST, on the other hand, end up worse than the HAR (N-BEATS with t = −5.7). Having more parameters gives a network no advantage if the structure it needs is that of a four-coefficient model.
• The iTransformer is the worst, for two reasons that can be named. Its advantage is relating series, and in vol5 each asset follows its own HAR (they share only a few shocks); and it sees six times fewer examples (one sample per date), so it overfits sooner. With real volatilities, correlated across assets, attention across series would have something to exploit.
• The library's N-HiTS falls behind the HAR: more parameters and an architecture designed for long horizons do not replace the right structure.

Exercises

1. The figures in the table come from one seed. Train N-BEATS (the fastest of the four to train) with seeds 0, 1 and 2 and report the mean and dispersion of the test error. Does the conclusion against the HAR change?

errores = []
for semilla in (0, 1, 2):
    torch.manual_seed(semilla)
    m = NBeats()
    entrenar(m, U, semilla=semilla)
    errores.append(np.mean((y_p - prever(m, U["prueba"][0])) ** 2))
print("MSE per seed:", np.round(errores, 3),
      f"mean {np.mean(errores):.3f}, std. {np.std(errores):.3f}")
print(f"HAR: {ecm_ref['HAR']:.3f}")

MSE per seed: [0.59  0.591 0.593] mean 0.591, std. 0.001
HAR: 0.577

2. DLinear and its relatives say that a linear model on the complete window is a tough rival. Fit by least squares a regression of the target on the 64 days of the window (65 coefficients) and compare it with the HAR (4 coefficients).

def con_cte(W):
    return np.c_[np.ones(len(W)), W.double().numpy()]

b64 = np.linalg.lstsq(con_cte(U["ent"][0]),
                      U["ent"][1].double().numpy(), rcond=None)[0]
e64 = np.mean((y_p - con_cte(U["prueba"][0]) @ b64) ** 2)
print(f"linear 64 days: {e64:.3f}   HAR: {ecm_ref['HAR']:.3f}")
print("weights of the last 5 days:", b64[-5:].round(3))
print("sum of the weights of days 23 to 64 back:",
      b64[1:43].sum().round(3))

linear 64 days: 0.580   HAR: 0.577
weights of the last 5 days: [0.082 0.086 0.061 0.096 0.32 ]
sum of the weights of days 23 to 64 back: 0.002

3. With real volatility data for 50 stocks of the same index, would you expect the iTransformer to improve its position relative to this lesson? What would you check before believing it?