Synthetic data for stress testing: GBM, GARCH and bootstrapping.

Why manufacture histories

A backtest on history answers one question: what would have happened with these data. Lesson 2.5.5 resampled that history to see how much the result could vary. But resampling only recombines what already happened: if the sample has no 35% drop in a month, no bootstrap is going to invent it. A synthetic data generator is a model that produces new series with whatever statistical properties are asked of it: it can be asked for a crisis that is not in the sample, a thousand years of data or a market with no edge at all.

Three legitimate uses and one that is not:

flowchart LR
  G[Generator] --> C[Test the code]
  G --> E[Stress tests]
  G --> R[Risk and sizing]
  G --> V[Validate an edge]
  C --> C2[Known answer]
  E --> E2[New crises]
  R --> R2[Tails and drawdown]
  V --> V2[No: circular]

• Test the code: if the generator plants an effect of known size, a correct method must recover it, and on a series with no effect it must find nothing. That is what all the data of this course do.
• Stress tests: insert scenarios that history does not have, or has only once, and see how the strategy responds.
• Risk and sizing: estimate the distribution of drawdowns or losses with more paths than history provides.
• Validate an edge: no. On synthetic data, a strategy wins if the generator contains what the strategy exploits. The reasoning is circular: we come back to this at the end.

Three generators

The reference series is the future on the level 2 fictitious index IDA, 26 years of sessions. It is itself synthetic, and that allows something that cannot be done with real data: comparing what each generator estimates with the model that actually produced the data (further below).

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

fut = pd.read_parquet("futuros2_continuos.parquet")
ida = (fut[fut["mercado"] == "IDA"].set_index("fecha")
       ["cierre_ajustado"])
r = np.log(ida).diff().dropna()
print(f"{len(r)} sessions, {r.index[0]:%Y-%m-%d} "
      f"to {r.index[-1]:%Y-%m-%d}")
print(f"Annual mean {252 * r.mean():.1%}, "
      f"volatility {np.sqrt(252) * r.std():.1%}")

6538 sessions, 2000-01-04 to 2025-12-31
Annual mean 6.3%, volatility 20.1%

Geometric Brownian motion

The simplest model, and the one behind Black-Scholes: the log of the price is a random walk with independent normal returns of constant mean $\mu$ and volatility $\sigma$:

N, T = 200, len(r)
rng = np.random.default_rng(2025)

mu, sigma = r.mean(), r.std()
gbm = mu + sigma * rng.standard_normal((N, T))
print(gbm.shape)

(200, 6538)

GJR-GARCH with t innovations

Bollerslev's GARCH (1986) makes today's variance depend on yesterday's shock and variance; the extension by Glosten, Jagannathan and Runkle (1993) adds a $\gamma$ term that only acts after negative days (the leverage effect), and Student's $t$ innovations give fatter tails than the normal (lesson 2.1.4):

arch fits the model by maximum likelihood. To simulate hundreds of paths it is faster to write the recursion with NumPy, one time step for all paths at once, than to call the model's simulate 200 times:

from arch import arch_model

ajuste = arch_model(100 * r, mean="Constant", vol="GARCH",
                    p=1, o=1, q=1, dist="t").fit(disp="off")
print(ajuste.params.round(3).to_dict())

def simular_gjr(p, T, N, rng, burn=500):
    """N paths of T sessions from a GJR-GARCH-t (as fractions)."""
    m, w, a, g, b, nu = (p["mu"], p["omega"], p["alpha[1]"],
                         p["gamma[1]"], p["beta[1]"], p["nu"])
    z = rng.standard_t(nu, (T + burn, N)) * np.sqrt((nu - 2) / nu)
    s2 = np.full(N, w / (1 - a - b - g / 2))
    e = np.zeros(N)
    out = np.empty((T + burn, N))
    for t in range(T + burn):
        s2 = w + (a + g * (e < 0)) * e ** 2 + b * s2
        e = np.sqrt(s2) * z[t]
        out[t] = m + e
    return out[burn:].T / 100

p_sim = ajuste.params.copy()
p_sim["mu"] = 100 * r.mean()     # the sample mean (see the text)
garch = simular_gjr(p_sim, T, N, rng)
persistencia = (ajuste.params["alpha[1]"] + ajuste.params["beta[1]"]
                + ajuste.params["gamma[1]"] / 2)
print(f"Persistence: {persistencia:.3f}; "
      f"mean of the fit: {252 * ajuste.params['mu'] / 100:.1%} annual")

{'mu': 0.054, 'omega': 0.032, 'alpha[1]': 0.044, 'gamma[1]': 0.104, 'beta[1]': 0.882, 'nu': 5.422}
Persistence: 0.978; mean of the fit: 13.5% annual

Two details. The persistence, $\alpha + \beta + \gamma/2$, is 0.978: a volatility shock takes about 31 sessions to halve ($\ln 0.5 / \ln 0.978$). And the estimated mean, 0.054% daily, is 13.5% annual, more than double the sample mean (6.3%). It is not a mistake: with $t$ innovations, the likelihood gives little weight to extreme days, which in IDA are mostly drops, and the estimated mean goes up. Simulating with it would draw a considerably more bullish market than the real one, so the simulation uses the sample mean.

Stationary bootstrap of the history

The third generator assumes no model: it reassembles history in blocks of random length, as in lesson 2.5.5, with an average of 75 sessions to preserve volatility streaks.

def indices_estacionario(n, largo, b, rng):
    idx = np.empty(largo, dtype=np.int64)
    idx[0] = rng.integers(n)
    salto = rng.random(largo) < 1 / b
    nuevo = rng.integers(n, size=largo)
    for t in range(1, largo):
        idx[t] = nuevo[t] if salto[t] else (idx[t - 1] + 1) % n
    return idx

x = r.to_numpy()
boot = np.array([x[indices_estacionario(T, T, 75, rng)]
                 for _ in range(N)])
print(gbm.shape, garch.shape, boot.shape)

(200, 6538) (200, 6538) (200, 6538)

What does each one reproduce?

The stylized facts of lesson 1.6.2, measured on IDA and on the median of each generator's 200 paths:

def acf(x, k):
    x = x - x.mean()
    return (x[k:] * x[:-k]).sum() / (x * x).sum()

def hechos(x):
    z = (x - x.mean()) / x.std()
    eq = np.exp(np.cumsum(x))
    return {"vol": np.sqrt(252) * x.std(),
            "curtosis": stats.kurtosis(x),
            "asimetría": stats.skew(x),
            "días >4σ": (np.abs(z) > 4).sum(),
            "acf r 1": acf(x, 1),
            "acf |r| 1": acf(np.abs(x), 1),
            "acf |r| 20": acf(np.abs(x), 20),
            "apalanc.": np.corrcoef(x[:-1], np.abs(x[1:]))[0, 1],
            "dd máx": (eq / np.maximum.accumulate(eq) - 1).min()}

tabla = {"IDA": hechos(x)}
for nombre, sims in (("GBM", gbm), ("GJR-t", garch),
                     ("bootstrap", boot)):
    tabla[nombre] = pd.DataFrame([hechos(s) for s in sims]).median()
pd.DataFrame(tabla).round(3)

               IDA    GBM   GJR-t  bootstrap
vol          0.201  0.201   0.182      0.201
curtosis    12.584 -0.004   9.009     12.141
asimetría   -0.373  0.001   0.105     -0.369
días >4σ    51.000  0.000  39.000     49.000
acf r 1      0.020  0.000  -0.003      0.019
acf |r| 1    0.247 -0.001   0.215      0.240
acf |r| 20   0.180  0.002   0.122      0.134
apalanc.    -0.056 -0.000  -0.038     -0.054
dd máx      -0.706 -0.540  -0.505     -0.661

GBM reproduces nothing but the mean and the volatility: no tails, no streaks, no skewness. GJR-GARCH-t reproduces short-term clustering (somewhat less: 0.21 against 0.25) and a good part of the tails and of the leverage effect (the $\gamma$ term provides it), but not the skewness of daily returns, because its $t$ innovations are symmetric. The bootstrap reproduces almost everything, as expected: they are the same days. Its median drawdown (−66.1%) comes close to the historical one (−70.6%); the other two stay at −54% and −51% because they spread the drops without clustering them.

The header figure shows the fact that none of them reproduces well: the long memory of volatility. IDA's autocorrelation of |r| is still 0.130 at 50 sessions and 0.077 at 100, and dies out around 150. That of GJR-t falls to 0.052 and 0.010, because its memory decays exponentially; that of the bootstrap, to 0.062 and 0.015, because the 75-session blocks cut any longer dependence.

RETARDOS = [1, 2, 5, 10, 20, 50, 100, 150, 200, 250]
memoria = pd.DataFrame(
    {"IDA": [acf(np.abs(x), k) for k in RETARDOS]}, index=RETARDOS)
for nombre, sims in (("GBM", gbm), ("GJR-t", garch),
                     ("bootstrap", boot)):
    memoria[nombre] = [np.median([acf(np.abs(s), k) for s in sims])
                       for k in RETARDOS]
memoria.rename_axis("retardo").round(3)

           IDA    GBM  GJR-t  bootstrap
retardo
1        0.247 -0.001  0.215      0.240
2        0.277  0.000  0.210      0.272
5        0.253  0.000  0.192      0.234
10       0.211  0.001  0.166      0.184
20       0.180  0.002  0.122      0.134
50       0.130  0.001  0.052      0.062
100      0.077 -0.001  0.010      0.015
150      0.007  0.002  0.002     -0.003
200      0.012 -0.001 -0.000     -0.002
250     -0.011  0.001 -0.005     -0.003

For a risk test, the consequence is direct: GBM and GJR-t place IDA's real drawdown in the tail of their distribution (below their 25th percentile: −60.6% in GBM and −60.2% in GJR-t), that is, they treat it as bad luck. The bootstrap leaves it inside its central box (25th percentile: −76.7%): it treats it as normal. Choosing the generator is choosing how much risk you are going to see.

What the fit does not give back

Since IDA is synthetic, it can be compared with the model that generated it. The level 2 README documents it: the equity market factor is a GJR-GARCH(1,1) with $t$ innovations ($\nu = 6$, $\alpha = 0.03$, $\gamma = 0.11$, $\beta = 0.86$ in the generator's code) whose long-run volatility and drift change according to a script of segments with fixed dates. From 2006, the IDA future follows exactly that factor (before, one of its own with the same script and the same parameters), plus the planted carry and momentum. The fit above returned $\alpha = 0.044$, $\gamma = 0.104$, $\beta = 0.882$ and $\nu = 5.4$: close in $\gamma$ and $\nu$, but with a persistence of 0.978 against the generator's 0.945 ($0.03 + 0.86 + 0.11/2$). The script can be read in the data, but not as is:

guion = {  # segments from the level 2 README: (from, to, vol, drift)
    "calma 2004-07": ("2004-01-02", "2007-07-13", 0.12, 0.07),
    "crisis 2008-09": ("2008-09-02", "2009-03-09", 0.52, -0.95),
    "calma 2012-15": ("2012-06-01", "2015-08-14", 0.12, 0.13),
    "covid 2020": ("2020-02-20", "2020-03-23", 0.80, -4.00),
    "2022": ("2022-01-04", "2022-10-12", 0.25, -0.32),
}
filas = {}
for k, (a, b, v, d) in guion.items():
    t = r.loc[a:b]
    filas[k] = {"sesiones": len(t), "vol_lp_guion": v,
                "vol_IDA": np.sqrt(252) * t.std(),
                "deriva_guion": d, "deriva_IDA": 252 * t.mean()}
pd.DataFrame(filas).T.round(2)

                sesiones  vol_lp_guion  vol_IDA  deriva_guion  deriva_IDA
calma 2004-07      888.0          0.12     0.13          0.07        0.09
crisis 2008-09     130.0          0.52     0.50         -0.95       -1.01
calma 2012-15      806.0          0.12     0.13          0.13        0.13
covid 2020          23.0          0.80     0.53         -4.00       -5.60
2022               195.0          0.25     0.19         -0.32       -0.05

In the long, quiet segments, realized volatility (0.13) comes close to the script's (0.12), and in the 130 sessions of 2008-2009 too (0.50 against 0.52). In the short 2020 segment, it falls short: the script sets the long-run volatility, and a GARCH takes dozens of sessions to approach it. In its 23 sessions, the script asked for 0.80 and realized volatility was 0.53. The drifts deviate more, because the GARCH-t noise sits on top of the script's: in 2020, −5.60 against −4.00; in 2022, −0.05 against −0.32 (the level 2 README documents, year by year, how far the realized factor deviates from the script). A GARCH fitted to the whole series does not see segments: it sees a volatility that rises and falls with a lot of inertia, and summarizes it in a higher persistence than the generator's (0.978 against 0.945).

It is a general lesson. The estimated parameters are those of the model being fitted, not those of the process that produced the data, and when the model is not the true one (with real data, it never is) the generator that comes out of the fit inherits its simplifications. To check that estimation does work when the model is the right one, it is enough to fit the GJR-t to its own paths:

import warnings
from arch.utility.exceptions import ConvergenceWarning

warnings.filterwarnings("ignore", category=ConvergenceWarning)
rng_ajuste = np.random.default_rng(7)
propios = simular_gjr(ajuste.params, T, 10, rng_ajuste)
reajustes = pd.DataFrame([
    arch_model(100 * s, mean="Constant", vol="GARCH", p=1, o=1, q=1,
               dist="t").fit(disp="off").params for s in propios])
pd.DataFrame({"verdad": ajuste.params,
              "media_10_ajustes": reajustes.mean(),
              "desv_10_ajustes": reajustes.std()}).round(3)

          verdad  media_10_ajustes  desv_10_ajustes
mu         0.054             0.059            0.007
omega      0.032             0.031            0.003
alpha[1]   0.044             0.047            0.009
gamma[1]   0.104             0.102            0.010
beta[1]    0.882             0.881            0.007
nu         5.422             5.351            0.351

On paths generated by the GJR-t itself, the fit recovers its parameters with a small error. On IDA it came close in $\alpha$, $\gamma$ and $\nu$, but not in persistence, because IDA is not a GJR-t with fixed parameters: its long-run volatility changes by segment. The ConvergenceWarning from arch (which the optimizer may emit on some series) is silenced by its category so that the output is stable; in real work it is worth reviewing.

Stress tests with inserted crises

The main use of a model-based generator is being able to insert scenarios. On the 200 GJR-t paths, starting at session 3,000, three crises of different shapes are inserted: a fast 35% drop in 23 sessions (like the 2020 one) with volatility tripled, a slow 45% drop over 500 sessions (like the 2000-2002 one) with volatility ×1.5, and a fast drop followed by a 45% rebound in 60 sessions (a V). In each segment the mean of the returns is replaced by the scenario's drift and their deviations are amplified:

T0 = 3000  # session at which the scenario starts

def escenario(sims, caida, sesiones, mult_vol, rebote=0.0,
              sesiones_rebote=0):
    """Insert a drop (and a rebound) into each simulated path."""
    out = sims.copy()
    tramo = slice(T0, T0 + sesiones)
    base = out[:, tramo] - out[:, tramo].mean(axis=1, keepdims=True)
    out[:, tramo] = np.log(1 + caida) / sesiones + mult_vol * base
    if sesiones_rebote:
        tramo2 = slice(T0 + sesiones, T0 + sesiones + sesiones_rebote)
        base2 = out[:, tramo2] - out[:, tramo2].mean(axis=1,
                                                     keepdims=True)
        out[:, tramo2] = (np.log(1 + rebote) / sesiones_rebote
                          + mult_vol * base2)
    return out

ESCENARIOS = {
    "base": garch,
    "caída rápida": escenario(garch, -0.35, 23, 3.0),
    "caída lenta": escenario(garch, -0.45, 500, 1.5),
    "en V": escenario(garch, -0.35, 23, 3.0, rebote=0.45,
                      sesiones_rebote=60),
}
for nombre, sims in ESCENARIOS.items():
    caida = np.exp(sims[:, T0:T0 + 560].cumsum(axis=1)).min(axis=1)
    print(f"{nombre:>13}: median worst point "
          f"{np.median(caida) - 1:.1%}")

         base: median worst point -10.1%
 caída rápida: median worst point -42.5%
  caída lenta: median worst point -50.0%
         en V: median worst point -37.2%

In the base scenario, the worst point of the 560 sessions after session 3,000 is, in median, 10.1% below the start. With the inserted crises, between −37.2% and −50.0%. Three simple strategies are put through each scenario: holding the index, holding it with a volatility target of 15% (position = 15% / 20-session EWMA volatility, capped at 2) and the 12-month trend of lesson 2.5.5 with the same scaling. All of them use only data up to the previous session.

def ewm_vol(sims, span=20):
    a = 2 / (span + 1)
    v = np.empty_like(sims)
    v[:, 0] = sims[:, :250].var(axis=1)
    for t in range(1, sims.shape[1]):
        v[:, t] = (1 - a) * v[:, t - 1] + a * sims[:, t - 1] ** 2
    return np.sqrt(252 * v)          # with data up to t - 1

def estrategias(sims):
    simple = np.expm1(sims)
    vol = ewm_vol(sims)
    tend = np.sign(pd.DataFrame(sims.T).rolling(252).sum()
                   .shift(1).fillna(0).to_numpy().T)
    peso_vol = np.minimum(0.15 / vol, 2.0)
    return {"comprar y mantener": simple,
            "objetivo de vol.": peso_vol * simple,
            "tendencia 12 m": tend * peso_vol * simple}

def peor_caida(rend):
    eq = np.cumprod(1 + rend[:, T0:T0 + 560], axis=1)
    return eq.min(axis=1) - 1

estres = {}
for esc, sims in ESCENARIOS.items():
    for est, rend in estrategias(sims).items():
        c = peor_caida(rend)
        estres[(esc, est)] = {"mediana": np.median(c),
                              "p5": np.quantile(c, 0.05)}
estres = pd.DataFrame(estres).T
for col in ("mediana", "p5"):
    print(f"Worst drawdown in 560 sessions ({col}):")
    print(estres[col].unstack().reindex(list(ESCENARIOS))
          .round(3).to_string(), end="\n\n")

Worst drawdown in 560 sessions (mediana):
              comprar y mantener  objetivo de vol.  tendencia 12 m
base                      -0.101            -0.085          -0.149
caída rápida              -0.425            -0.253          -0.185
caída lenta               -0.500            -0.401          -0.137
en V                      -0.372            -0.204          -0.228

Worst drawdown in 560 sessions (p5):
              comprar y mantener  objetivo de vol.  tendencia 12 m
base                      -0.442            -0.299          -0.398
caída rápida              -0.607            -0.433          -0.430
caída lenta               -0.651            -0.494          -0.265
en V                      -0.513            -0.336          -0.436

Each strategy has its bad scenario:

• Holding absorbs each crisis in full: −42.5% in the fast drop and −50.0% in the slow one (medians).
• The volatility target reduces the fast drop from −42.5% to −25.3%, because volatility spikes and the position shrinks within a few days; in the slow one it helps less (−40.1%), because volatility rises little and the position barely changes.
• Trend is the best in the slow drop (−13.7%): it has time to change sign and ends up short. In the V it loses more than in the fast drop with no rebound (−22.8% against −18.5%), because it goes short near the bottom and the rebound catches it the wrong way round. And in the base scenario, a GJR-t with no planted trend, it is the worst of the three (−14.9%): with no effect to capture, it only pays for the sign changes.

The 5th percentiles tell the same story more starkly: in the fast drop, one path in twenty leaves the holder with a loss of more than 60%. This is the information a stress test should give before trading: not whether the strategy wins, but in what kind of market it loses and how much.

The data of this course are synthetic

Everything above also describes how the data used in this course were made. The README of the level 2 data (linked in the resources) documents its generator, and it is a good example of what can and cannot be asked of a generator:

flowchart LR
  G[Segment script] --> M[GJR-t market]
  M --> A[Returns]
  S[GARCH-t sectors] --> A
  I[Idiosyncratic] --> A
  P[Planted effects] --> A
  E[Events] --> A
  A --> D[Files and README]

• A script of segments with fixed dates (calm, the 2000-2002 crisis, 2008, 2020, 2022…) sets the market's drift and long-run volatility; in the stress segments betas compress toward 1 and correlations rise.
• A GJR-GARCH market factor with $t$ innovations, GARCH-t sector factors and an idiosyncratic component per asset provide the tails, the clustering and the correlation structure.
• The planted effects have published formulas and sizes: the time-series momentum of the futures, the momentum, value and low-volatility premia of the stocks, the cointegration of the pairs, the carry of the currencies… Alongside them, controls with no effect: the “quality” and the ESG placebo of the stocks, the correlated but not cointegrated pairs, the ruptura series after its break, 995 of the 1,000 strategies in ruido2.
• The events (guided bankruptcies, splits, dividends, index rebalances, earnings) and the fundamentals are generated consistently with the prices.
• A fixed seed per dataset makes everything reproducible byte for byte with the command given in the README.

That construction is what makes it possible to use the data as a test bench: if a method applied to the futures momentum does not find it, the method fails; if it finds an effect in the ESG placebo more often than the significance level, it fails too. But it has limits that must be kept in mind every time a figure from this course is read:

The most important limit is the third one on the right. In real markets, published anomalies lose part of their return after publication11Note 1 · McLean, R. D. and Pontiff, J. (2016), “Does Academic Research Destroy Stock Return Predictability?”, The Journal of Finance 71(1), 5-32, doi:10.1111/jofi.12365: across 97 variables that predicted stock returns, portfolio returns are 26% lower out of the original sample and 58% lower after publication., and profitable strategies attract capital that erodes them. In the course data, the planted momentum is worth the same in 2001 as in 2025 because that is how it is written. No return in these data is evidence that a strategy works: it is evidence that the code finds what the generator contains.

Exercises

1. Add a fourth scenario: a 50% drop in 10 sessions with volatility multiplied by 4 (an extreme crash) and compute the median worst drawdown and the 5th percentile for the three strategies.

crash = escenario(garch, -0.50, 10, 4.0)
for est, rend in estrategias(crash).items():
    c = peor_caida(rend)
    print(f"{est:>18}: median {np.median(c):.1%}, "
          f"p5 {np.quantile(c, 0.05):.1%}")

comprar y mantener: median -55.5%, p5 -70.0%
  objetivo de vol.: median -27.0%, p5 -44.6%
    tendencia 12 m: median -18.2%, p5 -42.5%

2. Repeat the bootstrap with blocks of 250 sessions on average. Do the volatility memory at 100 sessions and the median drawdown change? What is gained and what is lost with such long blocks?

rng250 = np.random.default_rng(250)
boot250 = np.array([x[indices_estacionario(T, T, 250, rng250)]
                    for _ in range(N)])
h = pd.DataFrame([hechos(s) for s in boot250]).median()
print(f"ACF |r| at 100: "
      f"{np.median([acf(np.abs(s), 100) for s in boot250]):.3f}")
print(f"Median drawdown: {h['dd máx']:.1%}")

ACF |r| at 100: 0.041
Median drawdown: -63.7%

3. A colleague proposes validating an intraday reversal strategy by generating 10,000 days with a GJR-GARCH-t fitted to the index and checking that it wins in 90% of them. What would you tell them?