Lab: foundation models vs GARCH, HAR and LightGBM.

The test bench

Lesson 5.8.4 ended with a hint: the same Chronos-Bolt that does not predict returns explains almost half of the variance of the logarithm of realized variance. What remains is to know whether that is good, that is, whether it improves on what the usual volatility models already do. This lab measures it with five families of models:

The data are the synthetic vol5 data (level 5 data README, “Volatilidad” section). The true process is a HAR in the logarithm of integrated variance, with a leverage effect, jumps, level shifts and shocks on known dates. That has two consequences to keep in mind when reading the results: a well-estimated HAR is very close to optimal, and there is an oracle, the exact expectation of RV given the previous day's information, that sets the floor of what can be achieved. A model below the oracle would have a leak.

Three assets with different profiles are used: the index (strong leverage, γ = 0.18, no jumps), the bond (barely any leverage and a level shift of +0.60 in log variance on January 3, 2022) and the crypto asset (12 jumps a year with a 7% standard deviation and a level shift of −0.60 in 2019). The test covers 2018-2025.

flowchart TB
  D["vol5: 3 assets"] --> O["Origins 2018-2025"]
  O --> E["HAR and GARCH"]
  O --> G["LightGBM"]
  O --> F["Foundation models"]
  E --> P["QLIKE and MSE"]
  G --> P
  F --> P
  P --> T["DM and MCS"]

The protocol is a rolling origin without leakage. The models that are estimated (GARCH, HAR, LightGBM) are re-estimated every January 1 with all previous data and forecast every day of the year with those parameters; their features on day $t$ only use information up to $t-1$. The foundation models are not estimated: at each origin they receive the previous 512 sessions. Since each TimesFM forecast costs on the order of a tenth of a second on CPU, the models are compared on one session out of every 20 (101 origins per asset); at the end, what the cheap models say is checked with all sessions.

import numpy as np
import pandas as pd

d = pd.read_parquet("vol5_diario.parquet")
v = pd.read_parquet("vol5_verdad.parquet")
ACT = ["indice", "bono", "cripto"]
# overnight share ω, jumps per year and jump standard deviation (README)
PAR = {"indice": (0.20, 0, 0.0), "bono": (0.15, 0, 0.0),
       "cripto": (0.0, 12, 0.07)}
fechas = pd.DatetimeIndex(d.loc[d["activo"] == "indice", "fecha"])
n = len(fechas)
corte = fechas.searchsorted(pd.Timestamp("2018-01-02"))
orig = np.arange(corte, n, 20)              # one session out of every 20
inicio = {y: fechas.searchsorted(pd.Timestamp(f"{y}-01-01"))
          for y in range(2018, 2027)}
r = {a: d.loc[d["activo"] == a, "rendimiento"].to_numpy() for a in ACT}
rv = {a: d.loc[d["activo"] == a, "rv_5min"].to_numpy() for a in ACT}
lrv = {a: np.log(rv[a]) for a in ACT}
print(f"{len(orig)} origins per asset from {fechas[orig[0]].date()} "
      f"to {fechas[orig[-1]].date()}; {n - corte} sessions in 2018-2025")

101 origins per asset from 2018-01-02 to 2025-12-16; 2011 sessions in 2018-2025

Losses, oracle and HAR

The true variance is not observed; it is compared with the 5-minute RV, which measures it with error. Patton (2011) showed that not every loss works with a noisy proxy: MSE and QLIKE do rank the models the same as if the true variance were known. QLIKE, with $h$ the forecast,

is 0 when $h = RV$, penalizes falling short more than overshooting and depends only on the ratio $RV/h$, so a day of extreme variance does not dominate it the way it dominates MSE. The HAR is estimated in logarithms and its forecast of the mean carries the lognormal correction $e^{\hat x + \hat\sigma^2/2}$, with $\hat\sigma^2$ the variance of the training residuals. The oracle is the README formula: $E[RV_t \mid t-1] = (1-\omega)\,E[IV_t \mid t-1] + \lambda\,\sigma_J^2$, with the overnight share $\omega$, the jump intensity $\lambda$ (per session) and its variance $\sigma_J^2$.

flowchart LR
  A["January 1"] --> B["Re-estimate with < t"]
  B --> C["Forecast each day"]
  C --> D["Next year"]
  D --> A

def qlike(obs, f):
    """Patton's QLIKE: obs/f - log(obs/f) - 1 (0 if f = obs)."""
    return obs / f - np.log(obs / f) - 1

def medias_har(x):
    s = pd.Series(x)
    return pd.DataFrame({"d1": s.shift(1),
                         "s5": s.rolling(5).mean().shift(1),
                         "m22": s.rolling(22).mean().shift(1)})

prev, coef = {a: {} for a in ACT}, {}
for a in ACT:
    om, nj, sj = PAR[a]
    iv_prev = v.loc[v["activo"] == a, "iv_prevista"].to_numpy()
    prev[a]["oraculo"] = iv_prev * (1 - om) + nj / 252 * sj ** 2
    prev[a]["ingenua"] = np.r_[np.nan, rv[a][:-1]]
    X = medias_har(lrv[a]).to_numpy()
    har = np.full(n, np.nan)
    for y in range(2018, 2026):             # re-estimated every January
        t0, t1 = inicio[y], inicio[y + 1]
        ok = np.isfinite(X[:t0]).all(axis=1)
        A = np.column_stack([np.ones(ok.sum()), X[:t0][ok]])
        b, *_ = np.linalg.lstsq(A, lrv[a][:t0][ok], rcond=None)
        s2 = np.var(lrv[a][:t0][ok] - A @ b)
        har[t0:t1] = np.exp(b[0] + X[t0:t1] @ b[1:] + s2 / 2)
    prev[a]["har"], coef[a] = har, b
pd.DataFrame(coef, index=["c", "d1", "s5", "m22"]).round(3)

     indice   bono  cripto
c    -0.808 -1.215  -0.615
d1    0.373  0.273   0.137
s5    0.307  0.385   0.274
m22   0.238  0.236   0.503

The coefficients of the last fit (January 2025) show the HAR structure: the weight is spread among yesterday's RV, the weekly mean and the monthly mean, adding up to 0.92 on the index, 0.89 on the bond and 0.91 on the crypto asset, close to the true persistence (0.92). On the crypto asset the monthly mean weighs more (0.503): jumps put noise into the daily RV that the long mean filters out.

GARCH with arch

The GARCH(1,1) is estimated with daily returns, not with RV. arch_model(...).fit(last_obs=t0) estimates with the data before t0 and forecast(start=t0 - 1, ...) filters the rest of the year with those fixed parameters, so the forecast for day $t$ only uses returns up to $t-1$ (lesson 2.1.4). The GARCH variance is close to close; RV does not see the overnight part, so it is rescaled with the ratio $k$ between the mean RV and the variance of the training returns, as in the README.

from arch import arch_model

for a in ACT:
    gar = np.full(n, np.nan)
    for y in range(2018, 2026):
        t0, t1 = inicio[y], inicio[y + 1]
        am = arch_model(100 * r[a][:t1], mean="Constant", vol="GARCH",
                        p=1, q=1, dist="normal")
        res = am.fit(last_obs=t0, disp="off")
        # the row for t uses data up to t and forecasts t + 1
        h = res.forecast(start=t0 - 1, horizon=1, reindex=False)
        h = h.variance.to_numpy()[: t1 - t0, 0] / 1e4
        k = rv[a][:t0].mean() / np.var(res.resid[:t0] / 100)
        gar[t0:t1] = k * h        # to the RV scale (without the overnight part)
    prev[a]["garch"] = gar
print(res.params.round(4).to_dict(), f"k = {k:.2f}")

{'mu': 0.1113, 'omega': 0.1264, 'alpha[1]': 0.0183, 'beta[1]': 0.9706} k = 1.04

The output is from the last fit, the crypto asset in January 2025: α = 0.018 and β = 0.971, a persistence of 0.989 typical of a daily GARCH, and $k$ = 1.04 (the crypto asset has no overnight part: ω = 0). The GARCH reacts only to squared returns, a much noisier measure of variance than RV; that is why, on data with RV, the HAR usually beats it comfortably.

LightGBM with lags

The tree model receives more information than the HAR: five lags of log RV, the 5-, 22- and 66-session means, yesterday's return and the negative part of the returns of the last day and the last week, which is what it needs to capture the leverage effect (variance rises after declines). A small, deterministic configuration (lesson 5.5.3), re-estimated every January; the forecast of the mean carries the same lognormal correction as the HAR.

import lightgbm as lgb

def features(a):
    s, ret = pd.Series(lrv[a]), pd.Series(r[a])
    F = medias_har(lrv[a])
    for k in range(2, 6):
        F[f"d{k}"] = s.shift(k)
    F["m66"] = s.rolling(66).mean().shift(1)
    F["r1"] = ret.shift(1)
    F["neg1"] = np.minimum(ret, 0).shift(1)
    F["neg5"] = np.minimum(ret, 0).rolling(5).sum().shift(1)
    return F

def lgbm():
    return lgb.LGBMRegressor(
        n_estimators=300, learning_rate=0.03, num_leaves=8,
        min_child_samples=100, subsample=0.8, subsample_freq=1,
        colsample_bytree=0.8, random_state=0, n_jobs=4,
        deterministic=True, force_row_wise=True, verbose=-1)

ult = {}
for a in ACT:
    F = features(a)
    ok = F.notna().all(axis=1).to_numpy()
    p = np.full(n, np.nan)
    for y in range(2018, 2026):
        t0, t1 = inicio[y], inicio[y + 1]
        tr = ok & (np.arange(n) < t0)
        m = lgbm().fit(F[tr], lrv[a][tr])
        s2 = np.var(lrv[a][tr] - m.predict(F[tr]))
        p[t0:t1] = np.exp(m.predict(F[t0:t1]) + s2 / 2)
    prev[a]["lgbm"], ult[a] = p, m
imp = pd.Series(ult["indice"].booster_.feature_importance("gain"),
                index=F.columns)
(imp / imp.sum()).sort_values(ascending=False).head(5).round(3)

s5      0.437
d1      0.334
m22     0.096
r1      0.042
neg1    0.032
dtype: float64

The gain importance of the index model in 2025 is dominated by the same variables as the HAR: the weekly mean (43.7%), yesterday's RV (33.4%) and the monthly mean (9.6%). Yesterday's return (4.2%) and its negative part (3.2%) weigh little, but they are information a HAR-RV does not have; if LightGBM wins on some asset, that will be why.

The foundation models on log RV

Both models receive the logarithm of RV, not RV: in logarithms the series is almost Gaussian and its scale does not change with the level (exercise 1 tries RV in levels). Both return quantiles; the median of the log, exponentiated, is the median of RV, which sits below its mean because the distribution of RV has a right tail. To forecast the mean, which is what QLIKE and MSE reward, the lognormal correction is applied with the dispersion the model itself gives: if the log is approximately normal, $\sigma \approx (q_{0.9} - q_{0.1}) / (2 \cdot 1.2816)$ and the mean is $e^{q_{0.5} + \sigma^2/2}$. The 303 series (101 origins × 3 assets) go in a single batch.

import torch

torch.manual_seed(0)
torch.use_deterministic_algorithms(True)
from chronos import BaseChronosPipeline
import timesfm

CTX = 512
ctx = [lrv[a][t - CTX:t].astype("float32") for a in ACT for t in orig]
Z = 2 * 1.2815515655446004            # distance between q10 and q90 in σ

chronos = BaseChronosPipeline.from_pretrained(
    "amazon/chronos-bolt-small", device_map="cpu", dtype=torch.float32)
q, _ = chronos.predict_quantiles(
    inputs=[torch.tensor(c) for c in ctx], prediction_length=1,
    quantile_levels=[0.1, 0.5, 0.9])
q = q[:, 0, :].numpy().astype(float)
ch_med, ch_s = q[:, 1], (q[:, 2] - q[:, 0]) / Z

tfm = timesfm.TimesFM_2p5_200M_torch.from_pretrained(
    "google/timesfm-2.5-200m-pytorch", torch_compile=False)
tfm.compile(timesfm.ForecastConfig(
    max_context=CTX, max_horizon=16, normalize_inputs=True,
    use_continuous_quantile_head=True, per_core_batch_size=len(ctx)))
punto, cuant = tfm.forecast(horizon=1, inputs=ctx)
tf_med = punto[:, 0].astype(float)
tf_s = (cuant[:, 0, 9] - cuant[:, 0, 1]).astype(float) / Z

k = len(orig)
for i, a in enumerate(ACT):
    sl = slice(i * k, (i + 1) * k)
    for nom, med, s in [("chronos", ch_med, ch_s),
                        ("timesfm", tf_med, tf_s)]:
        prev[a][nom + "_mediana"] = np.exp(med[sl])
        prev[a][nom] = np.exp(med[sl] + s[sl] ** 2 / 2)
print("forecast log σ (mean):", round(ch_s.mean(), 3),
      round(tf_s.mean(), 3))

forecast log σ (mean): 0.409 0.42

The mean forecast dispersion of log RV is 0.41 for Chronos and 0.42 for TimesFM: the correction multiplies the median by $e^{0.41^2/2} \approx 1.09$, 9% more forecast variance.

Results: QLIKE and MSE

MOD = ["ingenua", "garch", "har", "lgbm", "chronos_mediana",
       "chronos", "timesfm_mediana", "timesfm", "oraculo"]

def en_orig(a, nom):
    f = prev[a][nom]
    return f if len(f) == len(orig) else f[orig]

perd = {a: {m: qlike(rv[a][orig], en_orig(a, m)) for m in MOD}
        for a in ACT}
tq = pd.DataFrame({a: {m: perd[a][m].mean() for m in MOD} for a in ACT})
ecm = pd.DataFrame({a: {m: np.mean((rv[a][orig] - en_orig(a, m)) ** 2)
                        for m in MOD} for a in ACT})
te = ecm / ecm.loc["oraculo"]
pd.concat({"QLIKE": tq, "ECM / ECM oráculo": te}, axis=1).round(3)

                 QLIKE               ECM / ECM oráculo
                indice   bono cripto            indice   bono cripto
ingenua          0.105  0.082  0.662             1.644  2.160  3.004
garch            0.135  0.119  0.362             1.790  1.843  1.013
har              0.080  0.056  0.383             1.191  0.995  1.029
lgbm             0.074  0.054  0.398             1.155  0.917  1.028
chronos_mediana  0.085  0.060  0.461             1.375  1.071  1.049
chronos          0.082  0.059  0.408             1.275  1.139  1.045
timesfm_mediana  0.089  0.060  0.465             1.324  1.025  1.049
timesfm          0.085  0.059  0.414             1.234  1.030  1.037
oraculo          0.069  0.054  0.350             1.000  1.000  1.000

The header figure summarizes the table by dividing each QLIKE by the HAR's on the same asset. What it says:

• The oracle is the floor (0.069, 0.054 and 0.350) and nobody crosses it in QLIKE: there is no leakage.
• The HAR sits close to the floor on the index and the bond, as the README anticipates: it is almost the true model.
• LightGBM improves on the HAR on the index (0.074 vs 0.080) and on the bond (0.054 vs 0.056), and falls behind on the crypto asset (0.398 vs 0.383).
• The two foundation models, without training anything, land between the naive forecast and the HAR: Chronos 0.082, 0.059 and 0.408; TimesFM 0.085, 0.059 and 0.414. They are 2-8% above the HAR and well below the naive forecast (0.105, 0.082 and 0.662). Chronos, with 48 million parameters, is slightly ahead of TimesFM, with 231 million.
• The mean correction matters, especially with jumps: on the crypto asset, Chronos's QLIKE drops from 0.461 with the median to 0.408 with the mean; TimesFM's, from 0.465 to 0.414.
• The GARCH is the worst on the index and the bond (0.135 and 0.119) and the best non-oracle model on the crypto asset (0.362). The section “A short sample misleads” comes back to this result.

MSE relative to the oracle shows why it is not used as the main loss: on the bond, LightGBM (0.917) and the HAR (0.995) come out below the oracle, something impossible in expectation. Over 101 origins MSE is decided by a few days of extreme variance; QLIKE, which looks at the ratio, is much more stable.

The figure shows the bond between 2021 and 2023. On January 3, 2022 its volatility rises permanently (+0.60 in log variance, ≈ 35% more volatility). Neither the HAR, re-estimated with the whole past, nor TimesFM, which normalizes each 512-session context, knows that the shift has happened; both absorb it through recent RV and follow almost the same path. The foundation model has no adaptation advantage over a HAR, which is also fed by the 1-, 5- and 22-day means.

Are the differences significant?

With 101 origins per asset, every difference has its error. First, Diebold-Mariano for each model against the HAR with QLIKE, with Newey-West variance (4 lags, although with origins 20 sessions apart the autocorrelation of the losses is small) and a two-sided test. A positive $t$ means the model loses more than the HAR.

from scipy import stats

def dm_hac(l1, l2, rezagos=4):
    """Diebold-Mariano t with Newey-West variance (two-sided)."""
    dd = l1 - l2
    u = dd - dd.mean()
    nn = len(dd)
    var = u @ u / nn
    for j in range(1, rezagos + 1):
        var += 2 * (1 - j / (rezagos + 1)) * (u[j:] @ u[:-j]) / nn
    t = dd.mean() / np.sqrt(var / nn)
    return t, 2 * stats.t.sf(abs(t), nn - 1)

COMP = ["ingenua", "garch", "lgbm", "chronos", "timesfm", "oraculo"]
pd.DataFrame({a: {m: "{:+.2f} ({:.2f})".format(
                     *dm_hac(perd[a][m], perd[a]["har"]))
                  for m in COMP} for a in ACT})

               indice          bono        cripto
ingenua  +2.33 (0.02)  +2.73 (0.01)  +1.93 (0.06)
garch    +3.74 (0.00)  +5.56 (0.00)  -0.19 (0.85)
lgbm     -0.97 (0.33)  -1.04 (0.30)  +0.83 (0.41)
chronos  +0.59 (0.56)  +1.33 (0.19)  +0.67 (0.50)
timesfm  +1.66 (0.10)  +1.14 (0.26)  +0.76 (0.45)
oraculo  -1.94 (0.06)  -1.01 (0.31)  -0.88 (0.38)

The foundation models lose to the HAR on all three assets, but no difference is significant: $t$ between +0.59 and +1.66 (TimesFM on the index, p = 0.10). LightGBM improves on the HAR on the index and the bond without significance ($t$ = −0.97 and −1.04). Only the GARCH on the index and the bond ($t$ = +3.74 and +5.56) and the naive forecast on the index and the bond are clearly worse. The oracle row is the most revealing: not even the perfect forecast beats the HAR significantly in this sample ($t$ = −1.94 on the index, p = 0.06). With 101 points, little can be claimed.

Diebold-Mariano compares pairs. To ask which models could be the best at the same time, without choosing a reference model, there is the Model Confidence Set of Hansen, Lunde and Nason (2011): it starts from all the models, tests whether they all have the same expected loss and, if that is rejected, removes the worst one and repeats. Those that remain form a set that contains the best model with probability $1-\alpha$. The version in the cell uses the $T_{max}$ statistic: for each model, its mean loss minus the mean of those still alive, divided by its standard error; the distribution under the null comes from a moving block bootstrap (blocks of 5, to respect the dependence). The arch library ships the same test (arch.bootstrap.MCS), which serves to check the hand-written implementation.

from arch.bootstrap import MCS

def mcs(L, alfa=0.10, B=2000, bloque=5, semilla=0):
    """Hansen, Lunde and Nason MCS with the T_max statistic and a
    moving block bootstrap. Returns the MCS p-value per model."""
    rng = np.random.default_rng(semilla)
    nn, _ = L.shape
    nb = -(-nn // bloque)
    ini = rng.integers(0, nn - bloque + 1, size=(B, nb))
    idx = (ini[:, :, None] + np.arange(bloque)).reshape(B, -1)[:, :nn]
    vivos, p_mcs, pmax = list(L.columns), {}, 0.0
    while len(vivos) > 1:
        X = L[vivos].to_numpy()
        di = X.mean(0) - X.mean()          # loss relative to the mean
        Xb = X[idx].mean(1)
        dib = Xb - Xb.mean(1, keepdims=True)
        se = np.sqrt(((dib - di) ** 2).mean(0))
        t = di / se
        Tb = ((dib - di) / se).max(1)       # T_max under the null
        pmax = max(pmax, float(np.mean(Tb >= t.max())))
        if pmax >= alfa:
            break
        peor = vivos[int(np.argmax(t))]
        p_mcs[peor] = pmax
        vivos.remove(peor)
    return {**p_mcs, **{m: 1.0 for m in vivos}}

CAND = ["ingenua", "garch", "har", "lgbm", "chronos", "timesfm"]
for a in ACT:
    L = pd.DataFrame({m: perd[a][m] for m in CAND})
    propio = mcs(L)
    ref = MCS(L, size=0.10, reps=2000, block_size=5, method="max",
              bootstrap="mbb", seed=0)
    ref.compute()
    dentro = sorted(m for m in CAND if propio[m] >= 0.1)
    print(f"{a:7s} own:  {dentro}\n"
          f"        arch: {sorted(ref.included)}")

indice  own:  ['chronos', 'har', 'lgbm']
        arch: ['chronos', 'har', 'lgbm']
bono    own:  ['chronos', 'har', 'lgbm', 'timesfm']
        arch: ['chronos', 'har', 'lgbm', 'timesfm']
cripto  own:  ['chronos', 'garch', 'har', 'ingenua', 'lgbm', 'timesfm']
        arch: ['chronos', 'garch', 'har', 'ingenua', 'lgbm', 'timesfm']

The hand-written implementation and the one in arch give the same sets at 10%. On the index, HAR, LightGBM and Chronos remain; on the bond, those three and TimesFM; on the crypto asset none can be ruled out, not even the naive forecast, whose QLIKE is 1.7 times the HAR's: with jumps of 7% standard deviation, a few days dominate the losses and 101 origins are not enough to separate anything. The MCS does not give a winner; it gives the models the data do not allow you to rule out.

A short sample misleads

The cheap models can be evaluated over the 2,011 sessions of 2018-2025 instead of over 101. The cell does so and repeats Diebold-Mariano for LightGBM against the HAR, now with 10 Newey-West lags because the sessions are consecutive.

BAR = ["ingenua", "garch", "har", "lgbm", "oraculo"]
todo = pd.DataFrame({a: {m: qlike(rv[a][corte:], prev[a][m][corte:])
                         .mean() for m in BAR} for a in ACT})
t_lg = {a: float(round(dm_hac(
            qlike(rv[a][corte:], prev[a]["lgbm"][corte:]),
            qlike(rv[a][corte:], prev[a]["har"][corte:]),
            rezagos=10)[0], 2)) for a in ACT}
print(f"LightGBM vs HAR t, {n - corte} sessions:", t_lg)
todo.round(3)

LightGBM vs HAR t, 2011 sessions: {'indice': -2.41, 'bono': 3.07, 'cripto': 1.33}
         indice   bono  cripto
ingenua   0.126  0.101   0.525
garch     0.149  0.113   0.385
har       0.089  0.067   0.366
lgbm      0.086  0.069   0.377
oraculo   0.080  0.065   0.325

Two rankings from the subsample flip. On the bond, LightGBM goes from 0.96 times the HAR's QLIKE to 1.03, and with 2,011 sessions the difference is significant against it ($t$ = +3.07). On the crypto asset, the GARCH goes from 0.95 to 1.05. On the index, on the other hand, LightGBM stays ahead (0.086 vs 0.089) and now with significance ($t$ = −2.41). The difference between the two assets has a plausible explanation: the index has a strong leverage effect (γ = 0.18) that LightGBM captures with the negative part of the returns and the HAR-RV does not see; on the bond (γ = 0.04) that information does not exist and the flexibility of the trees only adds variance. It is not a proof: the two return features weigh only 7% of the model's gain. The QLIKEs of the naive forecast, the HAR and the oracle match those in the README table (0.126, 0.089 and 0.080 on the index); the GARCH's differs a little because here it is re-estimated every year with arch.

The cost of each model

tam = {
    "HAR": 4, "GARCH(1,1)": 4,
    "LightGBM (hojas)": int(sum(
        t["num_leaves"] for t in
        ult["indice"].booster_.dump_model()["tree_info"])),
    "Chronos-Bolt small": sum(x.numel() for x in
                              chronos.model.parameters()),
    "TimesFM-2.5": sum(x.numel() for x in tfm.model.parameters()),
}
pd.Series(tam).map("{:,}".format)

HAR                             4
GARCH(1,1)                      4
LightGBM (hojas)            2,400
Chronos-Bolt small     47,718,016
TimesFM-2.5           231,289,280
dtype: str

The HAR and the GARCH have four parameters each; LightGBM, 2,400 leaves (300 trees of 8); Chronos-Bolt small, 47.7 million parameters, and TimesFM-2.5, 231.3 million. The figure plots size against the mean QLIKE relative to the HAR: fifty million times more parameters do not buy accuracy on a series whose process is a HAR.

Timing depends on the machine and its load, so the course does not run it; this is the code to measure it:

import time

t0 = time.perf_counter()
chronos.predict_quantiles(inputs=[torch.tensor(c) for c in ctx],
                          prediction_length=1, quantile_levels=[0.5])
print(f"Chronos, {len(ctx)} forecasts: "
      f"{time.perf_counter() - t0:.1f} s")
t0 = time.perf_counter()
tfm.forecast(horizon=1, inputs=ctx)
print(f"TimesFM, {len(ctx)} forecasts: "
      f"{time.perf_counter() - t0:.1f} s")

On the machine where the lesson was written (16 cores shared with other processes, load ≈ 12, 4 threads per process), the 24 HAR re-estimations (3 assets × 8 years) took 0.02 s; the GARCH ones with their forecasts, 0.9 s; the LightGBM ones, 1.9 s; the 303 Chronos-Bolt forecasts, 3.1 s, and the TimesFM-2.5 ones, 35 s (62 s with the machine more loaded). The underlying difference is not training but inference: the HAR forecasts with a dot product, and TimesFM needs a pass through 231 million parameters for each forecast, about 0.1-0.2 s on this CPU. Forecasting the 2,011 sessions of the three assets would take about 6,000 passes, on the order of 10-20 minutes here; on GPU (device="cuda") it would drop a lot, but it still would not compete with four multiplications.

What to do with this

It is the same conclusion as the Brini (2026) study with 50 real assets, cited in lesson 5.8.4: foundation models do not give a uniform improvement over a well-specified Log-HAR, and size or architecture matter less than having a serious econometric baseline.

Common mistakes

• Comparing with weak baselines. Against the naive forecast, the foundation models look excellent (22 to 38% less QLIKE); against the HAR, they lose.
• Using the median as the variance forecast. With QLIKE or MSE you forecast the mean; the median of a distribution with a right tail falls short, all the more with jumps.
• Trusting MSE with few observations. On the bond, two models “beat the oracle” in MSE.
• Drawing conclusions from 100 points. Not even the oracle is significantly better than the HAR; the ranking of LightGBM and the GARCH flips with the full sample.
• Leakage in the evaluation. Features that include day $t$, parameters estimated with the test year, normalizations with the whole series (lesson 5.2.4): the oracle sets a floor that no honest model crosses.
• Ignoring inference. A model with 231 million parameters costs at every forecast, not only when training.

Exercises

1. Repeat the Chronos-Bolt forecast feeding it RV in levels instead of its logarithm (and its median as the forecast). Which representation suits it better?

ctx_n = [rv[a][t - CTX:t].astype("float32") for a in ACT for t in orig]
qn, _ = chronos.predict_quantiles(
    inputs=[torch.tensor(c) for c in ctx_n], prediction_length=1,
    quantile_levels=[0.5])
qn = qn[:, 0, 0].numpy().astype(float)
pd.DataFrame({a: {"niveles": qlike(rv[a][orig],
                                   qn[i * k:(i + 1) * k]).mean(),
                  "log (mediana)": perd[a]["chronos_mediana"].mean(),
                  "log (media)": perd[a]["chronos"].mean()}
              for i, a in enumerate(ACT)}).round(3)

               indice   bono  cripto
niveles         0.096  0.061   0.499
log (mediana)   0.085  0.060   0.461
log (media)     0.082  0.059   0.408

2. The Mincer-Zarnowitz recalibration regresses the variable on the forecast: if the forecast is unbiased, the slope is 1 and the intercept 0. Compute, for each asset, the slope of the regression of log RV on the log of the HAR, LightGBM, Chronos and TimesFM forecasts at the origins. What does a slope below 1 indicate? Could that regression be used to correct the forecasts?

mz = {}
for a in ACT:
    yv = np.log(rv[a][orig])
    mz[a] = {m: np.polyfit(np.log(en_orig(a, m)), yv, 1)[0]
             for m in ["har", "lgbm", "chronos", "timesfm"]}
pd.DataFrame(mz).round(2)

         indice  bono  cripto
har        1.07  1.06    1.06
lgbm       1.11  1.11    1.06
chronos    1.01  0.96    0.92
timesfm    1.05  1.01    1.00

3. No code: your team uses a HAR to size positions with volatility targeting (lesson 4.2.2) across 40 futures. Someone proposes replacing it with TimesFM-2.5 “because it is the state of the art”. With what you have measured, what would you answer and what test would you ask for before deciding?