Estacionalidad de spreads: análisis multi-año, patrones y trampas.

Dos estacionalidades distintas

Que el crack de la gasolina «sube antes del verano» o que el gas de invierno «es caro» son frases con dos sentidos que conviene separar desde el principio.

• La forma estacional de la curva. Los contratos de entrega en verano de gasolina valen más que los de invierno; el gas de enero vale más que el de octubre. Es una diferencia de nivel entre meses de entrega, conocida por todos y descontada en los precios. La Energy Information Administration de EE. UU. explica las razones físicas: los márgenes de refino varían por estaciones y la demanda de gasolina suele subir en verano; el almacenamiento de gas cubre los picos de demanda de invierno (enlaces en los recursos). Comprar el contrato barato y vender el caro no gana esa diferencia: los dos contratos convergen cada uno a su propio precio de entrega.
• La deriva estacional. Lo que se puede operar es otra cosa: que el spread de dos contratos fijos tienda a subir en una parte del año y a bajar en otra, más de lo que la curva descontaba. Si existe, es la forma de ganar dinero con la estacionalidad.

Los datos del nivel separan las dos a propósito. Su modelo tiene una forma estacional de la curva (descontada) y, en cinco spreads, una deriva estacional plantada cuya suma acumulada es

con $d$ el día del año, $d^{\ast}$ el día del máximo y $A$ la amplitud en logaritmos: el spread sube del valle al máximo durante medio año y baja el otro medio. La curva no la descuenta, pero la atenúa en los contratos lejanos: un par de contratos con vencimiento a $\tau$ años la recibe multiplicada por $e^{-\kappa_Q \tau}$ ($\kappa_Q$ es la velocidad de reversión que descuenta la curva). Los otros cuatro spreads son controles: no tienen ninguna deriva estacional. Son datos sintéticos: ninguna cifra de esta lección dice que un spread real sea estacional ni cuánto.

import numpy as np
import pandas as pd
from scipy import stats

ver = pd.read_csv("spreads3_verdad.csv", index_col="spread")
NUEVE = ["crack_gasolina", "spark", "invierno_verano_gas", "crush",
         "trigo_maiz", "crack_gasoleo", "aceite_harina", "oro_plata",
         "calendario_crudo"]
v = ver.loc[NUEVE, ["patas", "control", "amplitud_estacional_log",
                    "fecha_pico", "kappa_q"]]
v.columns = ["patas", "control", "A", "pico", "kappa_q"]
v

                               patas  control      A   pico  kappa_q
spread
crack_gasolina               PEG/PEC    False  0.030  05-31      4.0
spark                        ELB/GNA    False  0.080  07-16      6.0
invierno_verano_gas  GNA ene/GNA oct    False  0.035  09-30      1.0
crush                    GRH/GRA/GRS    False  0.025  12-15      2.5
trigo_maiz                   GRT/GRM    False  0.018  01-14      1.0
crack_gasoleo                PEH/PEC     True    NaN    NaN      3.0
aceite_harina                GRA/GRH     True    NaN    NaN      0.8
oro_plata                    MTO/MTP     True    NaN    NaN      0.0
calendario_crudo         PEC 1.º/2.º     True    NaN    NaN      1.8

El crack de la gasolina tiene su máximo plantado el 31 de mayo; el spark (electricidad punta frente a gas), el 16 de julio; la prima de invierno del gas (enero frente a octubre), el 30 de septiembre; el crush de la soja, el 15 de diciembre, y el trigo frente al maíz, el 14 de enero. Los controles son el crack del gasóleo, el reparto del crush entre aceite y harina, el ratio oro/plata (un paseo aleatorio con un guion de episodios) y el calendario 1.º-2.º del crudo. Todos los spreads se miden en logaritmos: el de dos patas es $\ln F_1 - \ln F_2$, la apuesta de igual valor en las dos patas de la lección 3.2.4.

El análisis honesto tiene seis pasos, y cada uno evita una trampa de las que se ven más abajo:

flowchart TB
  A[Contratos fijos] --> B[Cambios diarios]
  B --> C[Armónico y F]
  C --> D[Calendario previo]
  D --> E[Sharpe y su error]
  E --> F[Pruebas múltiples]

Construir el spread sin saltos

La estacionalidad se mide sobre el cambio diario de un spread con los mismos contratos, nunca sobre una serie que cambia de contrato. Para un spread intercommodity (lección 3.2.2), las dos patas tienen que tener el mismo mes de entrega y cambiar de mes a la vez: se toma el primer mes común cuya fecha de roll no ha llegado en ninguna de las patas. El crush tiene tres patas (harina y aceite contra soja, con los pesos del margen de molienda), así que cada lado es una cesta.

cur = pd.read_parquet("futuros3_curvas.parquet")
ctr = pd.read_csv("futuros3_contratos.csv",
                  parse_dates=["vencimiento", "fecha_roll"])
px = cur.set_index(["fecha", "contrato"])["liquidacion"]
nombre = ctr.set_index(["mercado", "mes_entrega"])["contrato"]
fechas = pd.DatetimeIndex(np.sort(cur["fecha"].unique()))
f0 = np.r_[fechas[:1], fechas[:-1]]       # sesión anterior
LETRAS = "FGHJKMNQUVXZ"

def mes_comun(mercados, letras=None):
    """Mes de entrega que se tiene en todas las patas: el primero
    cuya fecha de roll no ha llegado en ninguna de las patas."""
    t = ctr[ctr["mercado"].isin(mercados)]
    if letras:
        mes = t["mes_entrega"].str[5:].astype(int)
        t = t[mes.map(lambda k: LETRAS[k - 1]).isin(list(letras))]
    g = t.groupby("mes_entrega").agg(n=("mercado", "nunique"),
                                     roll=("fecha_roll", "min"))
    roll = g.loc[g["n"] == len(mercados), "roll"].sort_index()
    roll = roll[::-1].cummin()[::-1]
    i = np.searchsorted(roll.to_numpy(), fechas.to_numpy(), "right")
    meses = np.r_[roll.index.to_numpy(), [None]]
    return pd.Series(meses[i], index=fechas)

def precio(f, contratos):
    return px.reindex(list(zip(f, contratos))).to_numpy()

def cesta(patas, f, mes):
    """Valor de una cesta [(mercado, peso)] con el mes `mes`."""
    v = np.zeros(len(f))
    for mk, peso in patas:
        c = nombre.reindex(list(zip([mk] * len(mes), mes))).to_numpy()
        v = v + peso * precio(f, c)
    return v

def simple(x):
    return x - 1

def cambio(largo, corto, letras=None, fn=np.log):
    """Cambio diario con los contratos del día anterior (sin saltos
    de roll): del log-spread ln(largo) - ln(corto) o, con fn=simple,
    rendimiento simple del largo menos el del corto."""
    mes = mes_comun({m for m, _ in largo + corto}, letras).shift(1)
    d = (fn(cesta(largo, fechas, mes) / cesta(largo, f0, mes))
         - fn(cesta(corto, fechas, mes) / cesta(corto, f0, mes)))
    return pd.Series(d, index=fechas).iloc[1:]

PATAS = {
    "crack_gasolina": ([("PEG", 42)], [("PEC", 1)], None),
    "crack_gasoleo": ([("PEH", 42)], [("PEC", 1)], None),
    "spark": ([("ELB", 1)], [("GNA", 7)], None),
    "crush": ([("GRH", 0.022), ("GRA", 0.11)], [("GRS", 0.01)],
              "FHKNQU"),
    "aceite_harina": ([("GRA", 1)], [("GRH", 1)], "FHKNQUVZ"),
    "oro_plata": ([("MTO", 1)], [("MTP", 1)], None),
    "trigo_maiz": ([("GRT", 1)], [("GRM", 1)], "HKNUZ"),
}
D = {k: cambio(*v) for k, v in PATAS.items()}               # log
R = {k: cambio(*v, fn=simple) for k, v in PATAS.items()}    # simple
print({k: int(v.notna().sum()) for k, v in D.items()})

{'crack_gasolina': 4023, 'crack_gasoleo': 4023, 'spark': 4023, 'crush': 4023, 'aceite_harina': 4023, 'oro_plata': 4023, 'trigo_maiz': 4023}

Las siete series tienen las 4.023 sesiones con cambio de 2010 a 2025. La regla del mes común importa: con la fecha de roll más tardía de las patas, en lugar de la más temprana, en los cracks una pata ya vencida deja unas veinte sesiones al año sin precio (el crudo deja de cotizar antes que la gasolina del mismo mes). La tabla de medidas del README de los datos usa esa otra regla, y por eso sus cifras de los cracks no coinciden con las de esta lección: la construcción del spread ya mueve el resultado.

Los dos calendarios se construyen aparte: el del crudo con el par cotizado de spreads3_cotizados.parquet (lección 3.2.5), y la prima de invierno del gas con el par fijo de enero frente al octubre anterior, que se mantiene hasta que vence octubre:

cot = pd.read_parquet("spreads3_cotizados.parquet")
pec = cot[cot["mercado"] == "PEC"].set_index("fecha")
c1 = pec["pata_cerca"].shift(1).reindex(fechas)
c2 = pec["pata_lejos"].shift(1).reindex(fechas)

def par(a, b, fn):
    return (fn(precio(fechas, a) / precio(f0, a))
            - fn(precio(fechas, b) / precio(f0, b)))

# Gas: enero del año siguiente frente a octubre, hasta que vence octubre
ano = np.where(fechas.month <= 10, fechas.year, fechas.year + 1)
ene = pd.Series([f"GNAF{(a + 1) % 100:02d}" for a in ano], fechas)
octb = pd.Series([f"GNAV{a % 100:02d}" for a in ano], fechas)
e1, o1 = ene.shift(1), octb.shift(1)
vivo = fechas <= ctr.set_index("contrato")["vencimiento"].reindex(o1)
for dic, fn in ((D, np.log), (R, simple)):
    dic["calendario_crudo"] = pd.Series(par(c1, c2, fn),
                                        index=fechas).iloc[1:]
    gas = np.where(vivo, par(e1, o1, fn), np.nan)
    dic["invierno_verano_gas"] = pd.Series(gas, index=fechas).iloc[1:]
D = {k: v.replace([np.inf, -np.inf], np.nan) for k, v in D.items()}
R = {k: v.replace([np.inf, -np.inf], np.nan) for k, v in R.items()}
pd.Series({k: v.std() * np.sqrt(252) for k, v in D.items()},
          name="vol_anual").round(3)

crack_gasolina         0.163
crack_gasoleo          0.140
spark                  0.287
crush                  0.103
aceite_harina          0.272
oro_plata              0.180
trigo_maiz             0.166
calendario_crudo       0.040
invierno_verano_gas    0.156
Name: vol_anual, dtype: float64

Las volatilidades van del 4 % anual del calendario del crudo al 29 % del spark. Importan para lo que viene: una deriva estacional de unos pocos puntos en un spread que se mueve un 20 % al año es difícil de ver con pocos años.

El gráfico estacional y su ruido

La herramienta clásica del análisis estacional es el perfil multi-año: el cambio acumulado del spread desde el 1 de enero, año a año, y su media. Es lo que publican los servicios de estacionalidad de futuros. Con los dos cracks:

def perfil(serie):
    """Cambio acumulado del log-spread desde el 1 de enero de cada
    año, medido cada semana: media y cuartiles entre años."""
    s = serie.fillna(0)
    acum = s.groupby(s.index.year).cumsum()
    semana = np.minimum((s.index.dayofyear - 1) // 7 + 1, 52)
    tabla = acum.groupby([s.index.year, semana]).last().unstack(0)
    return pd.DataFrame({"media": tabla.mean(axis=1),
                         "p25": tabla.quantile(0.25, axis=1),
                         "p75": tabla.quantile(0.75, axis=1),
                         "error": tabla.std(axis=1) / np.sqrt(16)})

PERF = {k: perfil(D[k]) for k in ["crack_gasolina", "crack_gasoleo"]}
for k, p in PERF.items():
    print(f"{k:15s} media: mínimo semana {p['media'].idxmin():2d} "
          f"({p['media'].min():+.3f}), máximo semana "
          f"{p['media'].idxmax():2d} ({p['media'].max():+.3f}); "
          f"error típico de la media en la semana 22: "
          f"{p.at[22, 'error']:.3f}")

crack_gasolina  media: mínimo semana 36 (-0.009), máximo semana 16 (+0.068); error típico de la media en la semana 22: 0.030
crack_gasoleo   media: mínimo semana 19 (-0.013), máximo semana 26 (+0.013); error típico de la media en la semana 22: 0.019

La media del crack de la gasolina sube un 6,8 % hasta la semana 16 (mediados de abril) y baja después; el máximo plantado está en la semana 22 (31 de mayo). El del gasóleo, sin estacionalidad, se queda entre −1,3 % y +1,3 %. El perfil apunta en la dirección correcta, pero dos cosas obligan a no fiarse de él a simple vista. La primera es su error: la media de 16 años tiene un error típico de 3 puntos en la semana 22, del mismo orden que la propia subida media en esas semanas. La segunda son los cuartiles de la figura: en los dos spreads, la mitad de los años quedan en una banda de unos ±10 puntos; los años individuales no siguen el perfil, y el control tiene una dispersión casi igual. Un gráfico de medias enseña una curva suave aunque los datos sean ruido, porque promediar suaviza; por eso hace falta un contraste.

Un armónico anual con su contraste

Un contraste natural para una deriva con forma de coseno es una regresión del cambio diario $\Delta s_t$ sobre los incrementos diarios de un coseno y un seno anuales:

La amplitud es $\sqrt{a^2 + b^2}$ y el máximo, el día $\operatorname{atan2}(b, a) \cdot 365{,}25 / 2\pi$. El contraste F de $a = b = 0$ pregunta si hay estacionalidad. Tiene dos ventajas sobre el gráfico: usa todos los días (no 16 años) y busca una forma concreta, la del modelo, en vez de cualquier curva.

W = 2 * np.pi / 365.25

def armonico(s):
    """Regresión del cambio diario sobre los incrementos diarios de
    cos y sin anuales: amplitud, día del máximo y p-valor F."""
    doy = s.index.dayofyear.to_numpy()
    dc = np.r_[np.nan, np.diff(np.cos(W * doy))]
    ds = np.r_[np.nan, np.diff(np.sin(W * doy))]
    ok = np.isfinite(s.to_numpy()) & np.isfinite(dc)
    X = np.c_[np.ones(ok.sum()), dc[ok], ds[ok]]
    y = s.to_numpy()[ok]
    b = np.linalg.lstsq(X, y, rcond=None)[0]
    rss1 = ((y - X @ b) ** 2).sum()
    rss0 = ((y - y.mean()) ** 2).sum()
    gl = len(y) - 3
    p = stats.f.sf((rss0 - rss1) / 2 / (rss1 / gl), 2, gl)
    pico = (np.arctan2(b[2], b[1]) / W) % 365.25
    return np.hypot(b[1], b[2]), pico, p

def fecha(dia):
    return (pd.Timestamp("2025-01-01")
            + pd.Timedelta(days=round(dia) - 1)).strftime("%m-%d")

filas = {}
for k in NUEVE:
    a, pico, p = armonico(D[k])
    filas[k] = {"control": ver.at[k, "control"],
                "A_plantada": ver.at[k, "amplitud_estacional_log"],
                "pico_plantado": ver.at[k, "fecha_pico"],
                "A_medida": a, "pico_medido": fecha(pico),
                "p_F": p}
med = pd.DataFrame.from_dict(filas, orient="index")
med.round(3)

                     control  A_plantada pico_plantado  A_medida pico_medido    p_F
crack_gasolina         False       0.030         05-31     0.032       05-21  0.000
spark                  False       0.080         07-16     0.052       08-17  0.001
invierno_verano_gas    False       0.035         09-30     0.016       10-16  0.095
crush                  False       0.025         12-15     0.013       12-28  0.042
trigo_maiz             False       0.018         01-14     0.021       01-12  0.033
crack_gasoleo           True         NaN           NaN     0.007       07-23  0.592
aceite_harina           True         NaN           NaN     0.010       07-01  0.767
oro_plata               True         NaN           NaN     0.016       06-11  0.204
calendario_crudo        True         NaN           NaN     0.001       09-02  0.915

Contra la respuesta conocida:

• El crack de la gasolina y el spark se recuperan. Máximos medidos el 21 de mayo y el 17 de agosto frente al 31 de mayo y el 16 de julio plantados, con p de 0,0002 y 0,0008. La amplitud medida del spark (0,052) es menor que la plantada (0,080): su $\kappa_Q = 6$ la atenúa en los contratos que se tienen. La del crack de la gasolina (0,032) es algo mayor que la plantada atenuada: el ruido también suma.
• El trigo-maíz y el crush se rozan. Máximos casi exactos (12 de enero frente a 14, 28 de diciembre frente a 15) y p de 0,033 y 0,042: significativos al 5 % uno a uno, por poco.
• La prima de invierno del gas no se ve. Máximo el 16 de octubre (plantado: 30 de septiembre), pero p = 0,095. Es un efecto real que el contraste no detecta con 16 años: un falso negativo.
• Los controles no dan nada. p entre 0,20 (oro-plata) y 0,91 (calendario del crudo), con amplitudes de 0,001 a 0,016. El oro-plata tiene una amplitud medida (0,016) mayor que la del crush (0,013), que sí tiene efecto: su volatilidad hace que el contraste no la tome en serio.

La trampa de la serie continua

La forma más rápida de hacer un perfil estacional es tomar la serie continua sin ajustar del primer contrato de cada pata (lección 1.4.6) y restar. Parece inofensivo. Con los mismos datos:

def frente(mk):
    """Precio del contrato que se tiene, sin ajustar (lección 1.4.6)."""
    cm = ctr[ctr["mercado"] == mk].sort_values("fecha_roll")
    i = np.searchsorted(cm["fecha_roll"], fechas, side="right")
    return pd.Series(precio(fechas, cm["contrato"].to_numpy()[i]),
                     index=fechas)

ing = {}
for k in ["crack_gasolina", "spark", "trigo_maiz",
          "crack_gasoleo", "aceite_harina"]:
    (a, _), (b, _) = PATAS[k][0][0], PATAS[k][1][0]
    nivel = np.log(frente(a)) - np.log(frente(b))
    A1, p1, f1 = armonico(nivel.diff())
    ing[k] = {"A_bien": med.at[k, "A_medida"],
              "pico_bien": med.at[k, "pico_medido"],
              "A_serie": A1, "pico_serie": fecha(p1), "p_serie": f1}
pd.DataFrame.from_dict(ing, orient="index").round(3)

                A_bien pico_bien  A_serie pico_serie  p_serie
crack_gasolina   0.032     05-21    0.090      04-30    0.000
spark            0.052     08-17    0.106      07-06    0.010
trigo_maiz       0.021     01-12    0.043      11-09    0.001
crack_gasoleo    0.007     07-23    0.004      07-18    0.911
aceite_harina    0.010     07-01    0.010      07-26    0.771

En los tres spreads con estacionalidad, la serie continua duplica o triplica la amplitud y mueve el máximo: el crack de la gasolina pasa de 0,032 a 0,090 con el máximo un mes antes (30 de abril), y el trigo-maíz, de enero a noviembre. Los contrastes salen más significativos que con los contratos fijos (p prácticamente 0 en el crack). La razón es la forma estacional de la curva: cada roll cambia de un mes de entrega a otro con distinto nivel estacional, y como los rolls caen en las mismas fechas cada año, sus saltos se repiten por calendario. El armónico los toma por una deriva. En los dos controles de la tabla los saltos del roll no tienen ese patrón anual y la serie continua no cambia el resultado (p de 0,91 y 0,77). Es el error más fácil de cometer y el más difícil de ver: la serie continua da un patrón más limpio y más significativo que el verdadero, y en otra fecha.

La misma confusión entre nivel y deriva aparece con el gas, donde la forma de la curva es enorme:

g = cur[cur["mercado"] == "GNA"].pivot(index="fecha",
                                         columns="contrato",
                                         values="liquidacion")
nivel, subida = [], []
for a in range(2010, 2025):
    ene_a, oct_a = f"GNAF{(a + 1) % 100:02d}", f"GNAV{a % 100:02d}"
    par = np.log(g[ene_a] / g[oct_a])
    nivel.append(par.dropna().mean())
    tramo = D["invierno_verano_gas"].loc[f"{a}-04-10":f"{a}-09-20"]
    subida.append(tramo.sum())
print(f"Nivel medio de ln(enero / octubre): {np.mean(nivel):+.3f}")
print(f"Subida media del 10-04 al 20-09: {np.mean(subida):+.3f} "
      f"(años positivos: {np.sum(np.array(subida) > 0)} de 15)")

Nivel medio de ln(enero / octubre): +0.170
Subida media del 10-04 al 20-09: +0.052 (años positivos: 12 de 15)

El contrato de enero vale de media un 17 % más que el octubre anterior (0,170 en log): es la prima de invierno que cualquier productor o almacenista conoce y que la curva descuenta. La subida media del spread entre abril y septiembre, el tramo que opera la estrategia, es de 0,052, con 12 años positivos de 15. Las dos cifras son estacionales, pero solo la segunda es una ganancia posible, y es un tercio de la primera.

Una estrategia por calendario y el error de su Sharpe

La estrategia más sencilla sobre una deriva estacional es mecánica: largo en el spread desde 10 días después del valle hasta 10 días antes del máximo, corto desde 10 días después del máximo hasta 10 días antes del valle, cada año, con la posición del día anterior. Los controles usan el calendario de la gasolina (máximo el 31 de mayo). El resultado de cada día es el de una posición con el mismo valor en las dos patas: el rendimiento simple de la pata larga menos el de la corta (la serie R de la celda de construcción); el armónico, en cambio, se estima con los cambios del log-spread (D). Con 16 años de datos hay, en la práctica, 16 observaciones de cada calendario, y el error típico de un Sharpe anual, con rendimientos aproximadamente normales e independientes (Lo, 2002; lección 2.5.4), es

con $T$ en años: unos 0,25 con 16 años, sea cual sea el spread.

def posicion(indice, pico):
    """+1 de valle + 10 días a máximo - 10; -1 de máximo + 10 a
    valle - 10; 0 el resto."""
    doy = indice.dayofyear.to_numpy()
    valle = (pico - 182.6) % 365.25
    sube = ((doy - (valle + 10)) % 365.25) < 162.6
    baja = ((doy - (pico + 10)) % 365.25) < 162.6
    return pd.Series(np.where(sube, 1.0, np.where(baja, -1.0, 0.0)),
                     index=indice)

def pico_de(k):   # los controles usan el calendario de la gasolina
    p = ver.at[k, "dia_pico"]
    return p if p == p else ver.at["crack_gasolina", "dia_pico"]

PNL = {}
filas = {}
for k in NUEVE:
    s = R[k]                    # rendimientos simples de las patas
    pnl = (posicion(s.index, pico_de(k)).shift(1) * s).fillna(0)
    PNL[k] = pnl
    anos = len(pnl) / 252
    sr = pnl.mean() / pnl.std() * np.sqrt(252)
    ee = np.sqrt((1 + sr**2 / 2) / anos)       # Lo (2002), anual
    por_ano = pnl.groupby(pnl.index.year).sum()
    filas[k] = {"control": ver.at[k, "control"], "sharpe": sr,
                "error_tipico": ee, "t": sr / ee,
                "anos_positivos":
                    f"{(por_ano > 0).sum()}/{len(por_ano)}"}
estr = pd.DataFrame.from_dict(filas, orient="index")
estr = estr.sort_values("sharpe", ascending=False)
estr.round(2)

                     control  sharpe  error_tipico     t anos_positivos
trigo_maiz             False    0.70          0.28  2.52          13/16
invierno_verano_gas    False    0.59          0.27  2.18          13/16
crack_gasolina         False    0.54          0.27  2.03          11/16
spark                  False    0.51          0.27  1.91          13/16
oro_plata               True    0.39          0.26  1.52           9/16
crush                  False    0.37          0.26  1.45          12/16
aceite_harina           True    0.32          0.26  1.25          11/16
crack_gasoleo           True    0.08          0.25  0.30           8/16
calendario_crudo        True   -0.08          0.25 -0.31           7/16

La figura de la cabecera dibuja la tabla con barras de dos errores típicos. Los cuatro primeros, trigo-maíz (0,70), prima de invierno del gas (0,59), crack de la gasolina (0,54) y spark (0,51), tienen estacionalidad plantada. El quinto es un control: el oro-plata, con 0,39, por encima del crush (0,37), que sí la tiene, y con solo 9 años positivos de 16. El otro control con calendario prestado, aceite-harina, da 0,32 con 11 años positivos de 16. Las barras se solapan casi todas: con 16 años, un Sharpe de 0,3 a 0,6 no distingue un efecto de una coincidencia. Fíjate también en lo que no se ve en el contraste F: la prima de invierno del gas, que el contraste no detectaba (p = 0,095), es aquí la segunda mejor; y el oro-plata, que el contraste descartaba (p = 0,20), parece prometedor. El Sharpe y el contraste miden cosas distintas: el F, si hay una forma anual concreta; el Sharpe, si una regla gana dinero, por la razón que sea.

Con nueve spreads probados a la vez, además, hay que corregir por pruebas múltiples (lección 2.1.8):

from statsmodels.stats.multitest import multipletests

t = estr["t"]
umbral = stats.norm.isf(0.05 / len(t))      # Bonferroni, unilateral
p_holm = multipletests(med["p_F"], alpha=0.05, method="holm")[1]
print(f"t que pide Bonferroni con 9 pruebas: {umbral:.2f}, "
      f"un Sharpe de {umbral * estr['error_tipico'].mean():.2f}")
pd.DataFrame({"control": estr["control"], "t": t.round(2),
              "pasa_bonferroni": t > umbral,
              "p_F": med["p_F"].reindex(estr.index).round(4),
              "p_F_holm": pd.Series(p_holm, med.index)
              .reindex(estr.index).round(3)})

t que pide Bonferroni con 9 pruebas: 2.54, un Sharpe de 0.67
                     control     t  pasa_bonferroni     p_F  p_F_holm
trigo_maiz             False  2.52            False  0.0327     0.229
invierno_verano_gas    False  2.18            False  0.0955     0.477
crack_gasolina         False  2.03            False  0.0002     0.002
spark                  False  1.91            False  0.0008     0.006
oro_plata               True  1.52            False  0.2039     0.816
crush                  False  1.45            False  0.0415     0.249
aceite_harina           True  1.25            False  0.7670     1.000
crack_gasoleo           True  0.30            False  0.5924     1.000
calendario_crudo        True -0.31            False  0.9145     1.000

Bonferroni con nueve pruebas exige una t de 2,54 (unilateral): un Sharpe de unos 0,67 con 16 años. Ninguno pasa; el trigo-maíz se queda en 2,52. Con la corrección de Holm sobre los p del contraste F, que usa los datos diarios y la forma del patrón, sí pasan dos: el crack de la gasolina (p corregido 0,002) y el spark (0,006). El resto, incluidos tres efectos reales, no se distinguen del azar con estos datos. Es la conclusión honesta de 16 años de datos con efectos de este tamaño, y es la que el análisis estacional de spreads suele saltarse.

Elegir el calendario sin mirar el futuro

En la tabla anterior, el calendario era el verdadero, que en datos reales nadie conoce. Lo normal es estimarlo: el día del máximo sale del perfil o del armónico. Si se estima con toda la muestra y se mide la estrategia en la misma muestra, el día del máximo es un parámetro ajustado a los datos (lección 2.5.1). La forma honesta es estimarlo cada año con los años anteriores, desde 2014 para tener al menos cuatro años:

def sharpe(x):
    return x.mean() / x.std() * np.sqrt(252)

fuera = {}
for k in NUEVE:
    s, lg = R[k], D[k]          # P&L simple; calendario con el log
    trozos = []
    for a in range(2014, 2026):
        _, pico, _ = armonico(lg[lg.index.year < a])
        sa = s[s.index.year == a]
        trozos.append(posicion(sa.index, pico).shift(1).fillna(0) * sa)
    wf = pd.concat(trozos).fillna(0)
    s2 = s[s.index.year >= 2014]
    _, pico_todo, _ = armonico(lg)
    dentro = (posicion(s2.index, pico_todo).shift(1) * s2).fillna(0)
    verdad = (posicion(s2.index, pico_de(k)).shift(1) * s2).fillna(0)
    fuera[k] = {"control": ver.at[k, "control"],
                "calendario_verdadero": sharpe(verdad),
                "estimado_con_todo": sharpe(dentro),
                "estimado_antes": sharpe(wf)}
pd.DataFrame.from_dict(fuera, orient="index").round(2)

                     control  calendario_verdadero  estimado_con_todo  estimado_antes
crack_gasolina         False                  0.30               0.39            0.24
spark                  False                  0.52               0.62            0.56
invierno_verano_gas    False                  0.51               0.47            0.68
crush                  False                  0.24               0.65            0.38
trigo_maiz             False                  0.94               0.97            0.83
crack_gasoleo           True                  0.32               0.29            0.06
aceite_harina           True                  0.13               0.38            0.05
oro_plata               True                  0.24               0.26            0.26
calendario_crudo        True                 -0.01              -0.20           -0.44

Las tres columnas cubren los mismos años, 2014-2025. Con el calendario estimado con toda la muestra, el crush pasa de 0,24 (calendario verdadero) a 0,65, y el control aceite-harina, de 0,13 a 0,38: la estimación busca las fechas que mejor encajan con el ruido de esos años. Con el calendario estimado solo con los años anteriores, el crush da 0,38 y el aceite-harina, 0,05; el crack del gasóleo, 0,06. La validación fuera de muestra quita la mayor parte de la ventaja falsa, pero no toda: el oro-plata da entre 0,24 y 0,26 en las tres columnas. No es un error del método. El ratio oro/plata de los datos sigue un guion de episodios (subidas y caídas largas con fecha fija) y, con 12 años, algunos caen en las mismas mitades del año. Fuera de muestra, un control puede ganar durante una década por coincidencia de calendario.

Lo que la estacionalidad exige en real

• Contratos fijos o series ajustadas. El perfil se calcula con el par que se tiene, sin saltos. Una serie continua sin ajustar fabrica estacionalidad en los mercados con forma estacional de curva, que son justo los estacionales.
• Pocas observaciones. Dieciséis años son dieciséis repeticiones de cada calendario. Con un Sharpe de 0,5 y 16 años, el intervalo al 95 % va de casi 0 a 1. Más años ayudan poco si el mercado ha cambiado (nuevas especificaciones, exportaciones de gas licuado, biocombustibles): el patrón de 1995 no tiene por qué ser el de 2025.
• Ventanas y márgenes. Los 10 días de margen en cada extremo son un parámetro más; cambiarlos cambia el resultado (ejercicio 1). Elegirlos mirando el resultado es otra forma de ajuste (lección 2.5.6).
• Costes. Una estrategia estacional entra y sale dos veces al año, y además rola: en un spread de calendario o intercommodity, cada operación paga las dos patas (lección 3.2.5).
• Riesgo de cola. Los spreads de energía y de electricidad tienen saltos (huracanes, olas de frío, paradas de refinería). El Sharpe de una estrategia estacional no enseña su peor año: hay que mirarlo (ejercicio 2).

Ejercicios

1. Repite la estrategia del crack de la gasolina y del control oro-plata con márgenes de 0, 10, 20 y 30 días en cada extremo de las dos mitades del año. ¿Cambia el orden entre los dos?

def posicion_margen(indice, pico, margen):
    doy = indice.dayofyear.to_numpy()
    valle = (pico - 182.6) % 365.25
    largo_tramo = 182.6 - 2 * margen
    sube = ((doy - (valle + margen)) % 365.25) < largo_tramo
    baja = ((doy - (pico + margen)) % 365.25) < largo_tramo
    return pd.Series(np.where(sube, 1.0, np.where(baja, -1.0, 0.0)),
                     index=indice)

tabla = {}
for k in ["crack_gasolina", "oro_plata"]:
    s = R[k]
    tabla[k] = {m: sharpe((posicion_margen(s.index, pico_de(k), m)
                           .shift(1) * s).fillna(0))
                for m in [0, 10, 20, 30]}
pd.DataFrame(tabla).T.round(2)

                  0     10    20    30
crack_gasolina  0.49  0.54  0.70  0.65
oro_plata       0.27  0.39  0.39  0.24

2. Para los nueve spreads, calcula el peor año de la estrategia (suma de sus rendimientos diarios) y el Sharpe en dos mitades, 2010-2017 y 2018-2025. ¿Cuántos mantienen el signo? ¿Hay diferencia entre spreads plantados y controles?

mitades = {}
for k in NUEVE:
    pnl = PNL[k]
    mitades[k] = {"control": ver.at[k, "control"],
                  "peor_ano": pnl.groupby(pnl.index.year).sum().min(),
                  "sharpe_2010_17": sharpe(pnl.loc[:"2017"]),
                  "sharpe_2018_25": sharpe(pnl.loc["2018":])}
pd.DataFrame.from_dict(mitades, orient="index").round(2)

                     control  peor_ano  sharpe_2010_17  sharpe_2018_25
crack_gasolina         False     -0.19            0.80            0.30
spark                  False     -0.53            0.67            0.34
invierno_verano_gas    False     -0.11            0.70            0.49
crush                  False     -0.13            0.21            0.53
trigo_maiz             False     -0.22            0.38            1.06
crack_gasoleo           True     -0.25           -0.22            0.38
aceite_harina           True     -0.49            0.46            0.17
oro_plata               True     -0.15            0.65            0.11
calendario_crudo        True     -0.08            0.19           -0.38

3. Diseña, sin código, el estudio de la estacionalidad del crack de la gasolina con datos reales de CME (gasolina RBOB y crudo WTI). ¿Qué datos pedirías, cómo construirías el spread, qué contrastarías y qué tendría que salir para operarlo?