Seasonality of spreads: multi-year analysis, patterns and traps.

Two different seasonalities

That the gasoline crack “rises before summer” or that winter gas “is expensive” are statements with two meanings that are worth separating from the start.

• The seasonal shape of the curve. Gasoline contracts for summer delivery are worth more than winter ones; January gas is worth more than October gas. It is a difference in level between delivery months, known to everyone and priced in. The U.S. Energy Information Administration explains the physical reasons: refining margins vary by season and gasoline demand usually rises in summer; gas storage covers the winter demand peaks (links in the resources). Buying the cheap contract and selling the expensive one does not earn that difference: each of the two contracts converges to its own delivery price.
• The seasonal drift. What can be traded is something else: that the spread of two fixed contracts tends to rise in one part of the year and fall in another, more than the curve priced in. If it exists, it is the way to make money from seasonality.

The data of the level separate the two on purpose. Its model has a seasonal shape of the curve (priced in) and, in five spreads, a planted seasonal drift whose cumulative sum is

with $d$ the day of the year, $d^{\ast}$ the day of the peak and $A$ the amplitude in logarithms: the spread rises from trough to peak for half a year and falls the other half. The curve does not price it in, but it dampens it in the far contracts: a pair of contracts expiring in $\tau$ years receives it multiplied by $e^{-\kappa_Q \tau}$ ($\kappa_Q$ is the mean-reversion speed priced in by the curve). The other four spreads are controls: they have no seasonal drift at all. This is synthetic data: no figure in this lesson says that a real spread is seasonal or by how much.

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

The gasoline crack has its planted peak on May 31; the spark spread (peak electricity versus gas), on July 16; the winter gas premium (January versus October), on September 30; the soybean crush, on December 15, and wheat versus corn, on January 14. The controls are the heating oil crack, the split of the crush between oil and meal, the gold/silver ratio (a random walk with a script of episodes) and the crude oil 1st-2nd calendar. All spreads are measured in logarithms: the two-leg one is $\ln F_1 - \ln F_2$, the equal-value bet on the two legs from lesson 3.2.4.

The honest analysis has six steps, and each one avoids one of the traps shown below:

flowchart TB
  A[Fixed contracts] --> B[Daily changes]
  B --> C[Harmonic and F]
  C --> D[Prior calendar]
  D --> E[Sharpe and its error]
  E --> F[Multiple testing]

Building the spread without jumps

Seasonality is measured on the daily change of a spread with the same contracts, never on a series that changes contract. For an intercommodity spread (lesson 3.2.2), the two legs must have the same delivery month and change month at the same time: you take the first common month whose roll date has not arrived in any of the legs. The crush has three legs (meal and oil against soybeans, with the weights of the crushing margin), so each side is a basket.

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]]       # previous session
LETRAS = "FGHJKMNQUVXZ"

def mes_comun(mercados, letras=None):
    """Delivery month held in all the legs: the first one whose
    roll date has not arrived in any of the legs."""
    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):
    """Value of a basket [(market, weight)] with month `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):
    """Daily change with the previous day's contracts (no roll
    jumps): of the log-spread ln(long) - ln(short) or, with fn=simple,
    simple return of the long side minus that of the short side."""
    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}

The seven series have all 4,023 sessions with a change from 2010 to 2025. The common-month rule matters: with the latest roll date of the legs, instead of the earliest, in the cracks an already expired leg leaves about twenty sessions a year without a price (crude oil stops trading before gasoline of the same month). The table of measures in the data's README uses that other rule, and that is why its crack figures do not match those of this lesson: the construction of the spread already moves the result.

The two calendars are built separately: the crude oil one with the listed pair from spreads3_cotizados.parquet (lesson 3.2.5), and the winter gas premium with the fixed pair of January versus the previous October, which is held until October expires:

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: January of the following year versus October, until October expires
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

Volatilities range from 4% a year for the crude oil calendar to 29% for the spark. They matter for what comes next: a seasonal drift of a few points in a spread that moves 20% a year is hard to see with few years.

The seasonal chart and its noise

The classic tool of seasonal analysis is the multi-year profile: the cumulative change of the spread since January 1, year by year, and its average. It is what futures seasonality services publish. With the two cracks:

def perfil(serie):
    """Cumulative change of the log-spread since January 1 of each
    year, measured every week: mean and quartiles across years."""
    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} mean: minimum week {p['media'].idxmin():2d} "
          f"({p['media'].min():+.3f}), maximum week "
          f"{p['media'].idxmax():2d} ({p['media'].max():+.3f}); "
          f"standard error of the mean in week 22: "
          f"{p.at[22, 'error']:.3f}")

crack_gasolina  mean: minimum week 36 (-0.009), maximum week 16 (+0.068); standard error of the mean in week 22: 0.030
crack_gasoleo   mean: minimum week 19 (-0.013), maximum week 26 (+0.013); standard error of the mean in week 22: 0.019

The average of the gasoline crack rises 6.8% until week 16 (mid-April) and falls afterwards; the planted peak is in week 22 (May 31). That of heating oil, with no seasonality, stays between −1.3% and +1.3%. The profile points in the right direction, but two things force you not to trust it at first sight. The first is its error: the 16-year average has a standard error of 3 points in week 22, of the same order as the average rise itself in those weeks. The second is the quartiles in the figure: in both spreads, half of the years fall within a band of about ±10 points; individual years do not follow the profile, and the control has almost the same dispersion. A chart of averages shows a smooth curve even if the data is noise, because averaging smooths; that is why a test is needed.

An annual harmonic with its test

A natural test for a cosine-shaped drift is a regression of the daily change $\Delta s_t$ on the daily increments of an annual cosine and sine:

The amplitude is $\sqrt{a^2 + b^2}$ and the peak, day $\operatorname{atan2}(b, a) \cdot 365.25 / 2\pi$. The F-test of $a = b = 0$ asks whether there is seasonality. It has two advantages over the chart: it uses every day (not 16 years) and it looks for a specific shape, the model's, instead of any curve.

W = 2 * np.pi / 365.25

def armonico(s):
    """Regression of the daily change on the daily increments of
    annual cos and sin: amplitude, day of the peak and F p-value."""
    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

Against the known answer:

• The gasoline crack and the spark are recovered. Measured peaks on May 21 and August 17 versus the planted May 31 and July 16, with p of 0.0002 and 0.0008. The measured amplitude of the spark (0.052) is smaller than the planted one (0.080): its $\kappa_Q = 6$ dampens it in the contracts held. That of the gasoline crack (0.032) is somewhat larger than the dampened planted one: noise adds up too.
• Wheat-corn and the crush come close. Almost exact peaks (January 12 versus 14, December 28 versus 15) and p of 0.033 and 0.042: significant at 5% one by one, barely.
• The winter gas premium is not seen. Peak on October 16 (planted: September 30), but p = 0.095. It is a real effect that the test does not detect with 16 years: a false negative.
• The controls give nothing. p between 0.20 (gold-silver) and 0.91 (crude oil calendar), with amplitudes from 0.001 to 0.016. Gold-silver has a measured amplitude (0.016) larger than that of the crush (0.013), which does have an effect: its volatility keeps the test from taking it seriously.

The continuous series trap

The fastest way to build a seasonal profile is to take the unadjusted continuous series of the front contract of each leg (lesson 1.4.6) and subtract. It looks harmless. With the same data:

def frente(mk):
    """Price of the contract held, unadjusted (lesson 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

In the three spreads with seasonality, the continuous series doubles or triples the amplitude and moves the peak: the gasoline crack goes from 0.032 to 0.090 with the peak a month earlier (April 30), and wheat-corn, from January to November. The tests come out more significant than with the fixed contracts (p practically 0 in the crack). The reason is the seasonal shape of the curve: each roll changes from one delivery month to another with a different seasonal level, and since rolls fall on the same dates every year, their jumps repeat with the calendar. The harmonic takes them for a drift. In the two controls in the table the roll jumps do not have that annual pattern and the continuous series does not change the result (p of 0.91 and 0.77). It is the easiest mistake to make and the hardest to see: the continuous series gives a cleaner and more significant pattern than the true one, and on another date.

The same confusion between level and drift appears with gas, where the shape of the curve is huge:

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"Mean level of ln(January / October): {np.mean(nivel):+.3f}")
print(f"Mean rise from 04-10 to 09-20: {np.mean(subida):+.3f} "
      f"(positive years: {np.sum(np.array(subida) > 0)} of 15)")

Mean level of ln(January / October): +0.170
Mean rise from 04-10 to 09-20: +0.052 (positive years: 12 of 15)

The January contract is worth on average 17% more than the previous October (0.170 in log): it is the winter premium that any producer or storage operator knows and that the curve prices in. The average rise of the spread between April and September, the segment the strategy trades, is 0.052, with 12 positive years out of 15. Both figures are seasonal, but only the second is a possible gain, and it is a third of the first.

A calendar strategy and the error of its Sharpe

The simplest strategy on a seasonal drift is mechanical: long the spread from 10 days after the trough until 10 days before the peak, short from 10 days after the peak until 10 days before the trough, every year, with the previous day's position. The controls use the gasoline calendar (peak on May 31). Each day's result is that of a position with the same value in both legs: the simple return of the long leg minus that of the short one (the R series from the construction cell); the harmonic, on the other hand, is estimated with the changes of the log-spread (D). With 16 years of data there are, in practice, 16 observations of each calendar, and the standard error of an annual Sharpe, with approximately normal and independent returns (Lo, 2002; lesson 2.5.4), is

with $T$ in years: about 0.25 with 16 years, whatever the spread.

def posicion(indice, pico):
    """+1 from trough + 10 days to peak - 10; -1 from peak + 10 to
    trough - 10; 0 otherwise."""
    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):   # the controls use the gasoline calendar
    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]                    # simple returns of the legs
    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), annual
    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

The header figure draws the table with two-standard-error bars. The first four, wheat-corn (0.70), the winter gas premium (0.59), the gasoline crack (0.54) and the spark (0.51), have planted seasonality. The fifth is a control: gold-silver, with 0.39, above the crush (0.37), which does have it, and with only 9 positive years out of 16. The other control with a borrowed calendar, oil-meal, gives 0.32 with 11 positive years out of 16. Almost all the bars overlap: with 16 years, a Sharpe of 0.3 to 0.6 does not distinguish an effect from a coincidence. Notice also what is not seen in the F-test: the winter gas premium, which the test did not detect (p = 0.095), is here the second best; and gold-silver, which the test ruled out (p = 0.20), looks promising. The Sharpe and the test measure different things: the F, whether there is a specific annual shape; the Sharpe, whether a rule makes money, for whatever reason.

With nine spreads tested at once, you also have to correct for multiple testing (lesson 2.1.8):

from statsmodels.stats.multitest import multipletests

t = estr["t"]
umbral = stats.norm.isf(0.05 / len(t))      # Bonferroni, one-sided
p_holm = multipletests(med["p_F"], alpha=0.05, method="holm")[1]
print(f"t required by Bonferroni with 9 tests: {umbral:.2f}, "
      f"a Sharpe of {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 required by Bonferroni with 9 tests: 2.54, a Sharpe of 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 with nine tests requires a t of 2.54 (one-sided): a Sharpe of about 0.67 with 16 years. None passes; wheat-corn stays at 2.52. With the Holm correction on the F-test p-values, which uses the daily data and the shape of the pattern, two do pass: the gasoline crack (corrected p 0.002) and the spark (0.006). The rest, including three real effects, cannot be distinguished from chance with this data. It is the honest conclusion from 16 years of data with effects of this size, and it is the one that seasonal spread analysis usually skips.

Choosing the calendar without looking at the future

In the previous table, the calendar was the true one, which nobody knows in real data. The usual thing is to estimate it: the day of the peak comes from the profile or from the harmonic. If it is estimated with the whole sample and the strategy is measured on the same sample, the day of the peak is a parameter fitted to the data (lesson 2.5.1). The honest way is to estimate it every year with the previous years, from 2014 to have at least four years:

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

fuera = {}
for k in NUEVE:
    s, lg = R[k], D[k]          # simple P&L; calendar with the 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

The three columns cover the same years, 2014-2025. With the calendar estimated on the whole sample, the crush goes from 0.24 (true calendar) to 0.65, and the oil-meal control, from 0.13 to 0.38: the estimation looks for the dates that best fit the noise of those years. With the calendar estimated only with the previous years, the crush gives 0.38 and oil-meal, 0.05; the heating oil crack, 0.06. Out-of-sample validation removes most of the false edge, but not all of it: gold-silver gives between 0.24 and 0.26 in all three columns. It is not a flaw of the method. The gold/silver ratio in the data follows a script of episodes (long rises and falls on fixed dates) and, with 12 years, some fall in the same halves of the year. Out of sample, a control can win for a decade by calendar coincidence.

What seasonality demands in live trading

• Fixed contracts or adjusted series. The profile is computed with the pair held, without jumps. An unadjusted continuous series manufactures seasonality in markets with a seasonal curve shape, which are precisely the seasonal ones.
• Few observations. Sixteen years are sixteen repetitions of each calendar. With a Sharpe of 0.5 and 16 years, the 95% interval goes from almost 0 to 1. More years help little if the market has changed (new specifications, liquefied natural gas exports, biofuels): the 1995 pattern need not be the 2025 one.
• Windows and buffers. The 10-day buffers at each end are one more parameter; changing them changes the result (exercise 1). Choosing them by looking at the result is another form of fitting (lesson 2.5.6).
• Costs. A seasonal strategy enters and exits twice a year, and it also rolls: in a calendar or intercommodity spread, each trade pays for both legs (lesson 3.2.5).
• Tail risk. Energy and electricity spreads have jumps (hurricanes, cold snaps, refinery outages). The Sharpe of a seasonal strategy does not show its worst year: you have to look at it (exercise 2).

Exercises

1. Repeat the gasoline crack strategy and the gold-silver control with buffers of 0, 10, 20 and 30 days at each end of the two halves of the year. Does the ranking between the two change?

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. For the nine spreads, compute the strategy's worst year (sum of its daily returns) and the Sharpe in two halves, 2010-2017 and 2018-2025. How many keep their sign? Is there a difference between planted spreads and controls?

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. Design, without code, the study of the seasonality of the gasoline crack with real CME data (RBOB gasoline and WTI crude). What data would you request, how would you build the spread, what would you test and what would have to come out to trade it?