Volatility models: ARCH, GARCH, EGARCH and GJR with the arch library.

From the EWMA to the GARCH

Lesson 1.6.1 estimated volatility with an exponential moving average (EWMA) of squared returns, and 1.6.2 showed why it works: the size of moves has memory (volatility clustering) and rises more after drops (leverage effect). Models of the GARCH family formalize those two ideas and add a third one that the EWMA does not have: volatility returns to a long-run level.

With $\varepsilon_t = r_t - \mu$ the return minus its mean and $\sigma_t^2$ its conditional variance (the one forecast for $t$ with the information up to $t-1$):

• $\alpha$: how much yesterday's news ($\varepsilon_{t-1}^2$) moves today's variance.
• $\beta$: how much of yesterday's variance carries over.
• $\alpha + \beta < 1$: persistence. The variance returns to its long-run level $\sigma^2_{LP} = \omega/(1-\alpha-\beta)$, and the deviation from it is multiplied by $\alpha + \beta$ every day.

The EWMA with parameter $\lambda$ is the case $\omega = 0$, $\alpha = 1-\lambda$, $\beta = \lambda$: persistence 1, with no level to return to. The variants in this lesson change the reaction to news (symmetric or not) and the distribution of the shocks:

flowchart TB
  A[ARCH q] --> G[GARCH 1,1]
  G --> J[GJR: asymmetric]
  G --> E[EGARCH: log var]
  G --> W[EWMA: no level]
  G --> D[Shocks N or t]

• ARCH(q) (Engle, 1982)11Note 1 · R. F. Engle, “Autoregressive Conditional Heteroscedasticity with Estimates of the Variance of United Kingdom Inflation”, Econometrica 50(4), 1982, 987-1007, https://doi.org/10.2307/1912773; D. B. Nelson, “Conditional Heteroskedasticity in Asset Returns: A New Approach”, Econometrica 59(2), 1991, 347-370, https://doi.org/10.2307/2938260; R. F. Engle and V. K. Ng, “Measuring and Testing the Impact of News on Volatility”, Journal of Finance 48(5), 1993, 1749-1778, https://doi.org/10.1111/j.1540-6261.1993.tb05127.x.: only the last $q$ values of $\varepsilon^2$. It needs many lags for that much memory; GARCH (Bollerslev, 1986) achieves it with one.
• GJR-GARCH (Glosten, Jagannathan and Runkle, 1993): adds $\gamma\,I_{t-1}\,\varepsilon_{t-1}^2$, with $I_{t-1} = 1$ if $\varepsilon_{t-1} < 0$. Drops weigh $\alpha + \gamma$; rises, $\alpha$.
• EGARCH (Nelson, 1991): models $\ln\sigma_t^2$, so it needs no sign restrictions. In the arch version, with $e_t = \varepsilon_t/\sigma_t$:

A negative $\gamma$ makes drops raise volatility more.

An index with known volatility

The data are synthetic (level 2 README). An equal-weighted index is built with the 120 fictional stocks of universo2_ajustados.parquet: the daily mean of the total returns of those that trade in each session, from 2006 to 2025. Its volatility is dominated by the generator's market factor, a GJR-GARCH(1,1) with Student's t innovations with 6 degrees of freedom and parameters $\alpha = 0.03$, $\gamma = 0.11$, $\beta = 0.86$ (curso-code/make_data_nivel2.py), whose long-run volatility changes in stretches: from 12% a year in calm periods to 52% in the 2008 crisis and 80% in March 2020 (README table). On top of that come sector and idiosyncratic factors with their own GARCH, which diversification dampens. The known answer is, therefore, approximate: that of a factor that weighs a lot, not that of the whole series.

import numpy as np
import pandas as pd
from statsmodels.stats.diagnostic import het_arch

aj = pd.read_parquet("universo2_ajustados.parquet")
px = aj.pivot(index="fecha", columns="ticker",
              values="cierre_ajustado")
r = 100 * px.pct_change(fill_method=None).mean(axis=1).dropna()

print(f"{len(r)} sessions, from {r.index[0]:%Y-%m-%d} "
      f"to {r.index[-1]:%Y-%m-%d}")
print(f"annual volatility: {r.std() * np.sqrt(252):.1f}%; "
      f"excess kurtosis: {r.kurt():.1f}")
lm = het_arch(r - r.mean(), nlags=5, result_object=True)
print(f"ACF of r² at lag 1: {(r**2).autocorr(1):.3f}; "
      f"ARCH-LM test (5 lags): {lm.lm:.0f}, "
      f"p = {lm.lmpval:.1e}")

5030 sessions, from 2006-01-04 to 2025-12-31
annual volatility: 20.1%; excess kurtosis: 10.8
ACF of r² at lag 1: 0.160; ARCH-LM test (5 lags): 688, p = 2.1e-146

Annual volatility of 20% and a kurtosis of 10.8: fat tails. The autocorrelation of $r^2$ at lag 1 is 0.16 and Engle's ARCH-LM test (regression of $\varepsilon_t^2$ on its five lags) rejects with a p-value on the order of $10^{-146}$: today's variance depends on yesterday's. result_object=True asks for a result with names (lm, lmpval); without it, statsmodels warns with a FutureWarning that the current tuple will change. It is the starting point of any GARCH.

Four models with arch

arch_model separates the mean (mean="Constant"), the volatility (vol="GARCH" with p lags of $\varepsilon^2$, o asymmetric ones and q of $\sigma^2$; or vol="EGARCH") and the distribution (dist="normal" or "t"). It is estimated by maximum likelihood. It is best to work in percent: with returns as decimals, $\omega$ is on the order of $10^{-6}$ and the optimizer struggles (the library warns).

from arch import arch_model

especs = {
    "GARCH-N": dict(vol="GARCH", p=1, q=1, dist="normal"),
    "GARCH-t": dict(vol="GARCH", p=1, q=1, dist="t"),
    "GJR-t": dict(vol="GARCH", p=1, o=1, q=1, dist="t"),
    "EGARCH-t": dict(vol="EGARCH", p=1, o=1, q=1, dist="t"),
}
ajustes = {k: arch_model(r, mean="Constant", **v).fit(disp="off")
           for k, v in especs.items()}
tabla = pd.DataFrame({k: f.params for k, f in ajustes.items()}).T
tabla["logL"] = [f.loglikelihood for f in ajustes.values()]
tabla["AIC"] = [f.aic for f in ajustes.values()]
tabla["BIC"] = [f.bic for f in ajustes.values()]
cols = ["omega", "alpha[1]", "gamma[1]", "beta[1]", "nu", "logL",
        "BIC"]
tabla[cols].round(3)

          omega  alpha[1]  gamma[1]  beta[1]     nu      logL        BIC
GARCH-N   0.039     0.100       NaN    0.871    NaN -7322.123  14678.339
GARCH-t   0.035     0.094       NaN    0.880  5.622 -7137.193  14317.001
GJR-t     0.032     0.033     0.100    0.892  5.888 -7110.901  14272.941
EGARCH-t  0.010     0.174    -0.068    0.979  5.759 -7122.623  14296.385

Three readings of the table:

• Student's t matters more than anything else. Going from normal to t innovations improves the log-likelihood by 185 points with just one more parameter ($\nu = 5.6$ degrees of freedom). Even with changing volatility, the standardized shocks have fat tails: the conclusion of lesson 1.6.2, now inside the model.
• The asymmetry is clear. The GJR improves on the GARCH-t by 26 log-likelihood points; its $\gamma = 0.100$ is close to the market factor's 0.11, and its $\alpha = 0.033$, to the 0.03. The EGARCH, with $\gamma = -0.068$, also captures it, although it fits worse than the GJR (BIC of 14,296 versus 14,273).
• $\beta$ and the persistence come out higher than in the generator. This is examined in the next section.

A likelihood ratio test formalizes the asymmetry: the GARCH-t is the GJR-t with $\gamma = 0$, so $2(\ln L_{GJR} - \ln L_{GARCH})$ is distributed as a $\chi^2$ with one degree of freedom if the null is true.

from scipy import stats

lr = 2 * (ajustes["GJR-t"].loglikelihood
          - ajustes["GARCH-t"].loglikelihood)
print(f"LR = {lr:.1f}; p = {stats.chi2.sf(lr, 1):.1e}")

LR = 52.6; p = 4.1e-13

Persistence and half-life

If today the variance is above its long-run level, the $h$-day forecast is

with $P$ the persistence: $\alpha + \beta$ in the GARCH, $\alpha + \beta + \gamma/2$ in the GJR (with symmetric shocks, half are negative) and $\beta$ in the EGARCH (for the log of the variance). The half-life of a shock is the $h$ that leaves the deviation at half: $\ln 0.5 / \ln P$.

def persistencia(f):
    p = f.params
    if "EGARCH" in f.model.volatility.name:
        return p["beta[1]"]
    return p["alpha[1]"] + p["beta[1]"] + p.get("gamma[1]", 0) / 2

filas = {}
for k, f in ajustes.items():
    P = persistencia(f)
    lp = (np.nan if "EGARCH" in k
          else np.sqrt(252 * f.params["omega"] / (1 - P)))
    filas[k] = {"persistencia": P, "vol_LP_%": lp,
                "vida_media": np.log(0.5) / np.log(P)}
filas["market (generator)"] = {
    "persistencia": 0.03 + 0.86 + 0.11 / 2, "vol_LP_%": np.nan,
    "vida_media": np.log(0.5) / np.log(0.945)}
pd.DataFrame(filas).T.round(3)

                    persistencia  vol_LP_%  vida_media
GARCH-N                    0.971    18.420      23.289
GARCH-t                    0.974    18.404      26.563
GJR-t                      0.976    18.224      28.090
EGARCH-t                   0.979       NaN      32.829
market (generator)         0.945       NaN      12.253

The estimated GJR has a persistence of 0.976 and a half-life of 28 sessions; the generator's market factor, 0.945 and 12 sessions. The model is not wrong because of poor estimation: the series has two sources of persistence and the GARCH has only one parameter for both. One is the dynamics of the shocks (the generator's). The other is the regimes: the long-run volatility of the script changes in stretches and holds for months or years, and a GARCH with a single $\sigma^2_{LP}$ interprets it as shocks that take a long time to dissipate. It is the result of Lamoureux and Lastrapes (1990): structural changes in the unconditional variance inflate the estimated persistence.

It is checked by estimating only on stretches where the script keeps the long-run volatility constant (12% a year):

tramos = {"2012-06 to 2015-08": ("2012-06-01", "2015-08-14"),
          "2016-02 to 2018-09": ("2016-02-12", "2018-09-28"),
          "2006-2025 full": ("2006-01-01", "2025-12-31")}
filas = {}
for k, (a, b) in tramos.items():
    f = arch_model(r[a:b], mean="Constant", **especs["GJR-t"]).fit(
        disp="off")
    P = persistencia(f)
    filas[k] = {"n": len(r[a:b]), "alpha": f.params["alpha[1]"],
                "gamma": f.params["gamma[1]"],
                "beta": f.params["beta[1]"], "persistencia": P,
                "vida_media": np.log(0.5) / np.log(P)}
pd.DataFrame(filas).T.round(3)

                         n  alpha  gamma   beta  persistencia  vida_media
2012-06 to 2015-08   806.0  0.009  0.069  0.920         0.964      18.964
2016-02 to 2018-09   663.0  0.000  0.153  0.882         0.958      16.261
2006-2025 full      5030.0  0.033  0.100  0.892         0.976      28.090

In the two calm stretches, the persistence drops to 0.964 and 0.958 (half-lives of 19 and 16 sessions), closer to the generator's 0.945. With 663 and 806 sessions, the individual parameters are noisy ($\alpha$ reaches zero in one of them, and $\gamma$ ranges from 0.07 to 0.15), but the direction is clear. With real data the stretches are not known; that is why there are models with a slow volatility component (for example, two-component models or regime-switching models, lesson 2.1.7).

The next figure compares the conditional volatility of the GJR, subsampled every 10 sessions, with the long-run volatility of the market factor in each stretch of the script:

News impact curve

The news impact curve (Engle and Ng, 1993) answers a concrete question: if today the variance is at its long-run level and tomorrow the index moves $x$%, what volatility does the model forecast for the day after tomorrow? It is the most direct way to see the asymmetry, and it is the figure at the top.

def impacto(f, eps):
    p = f.params
    if "EGARCH" in f.model.volatility.name:
        s2 = np.exp(p["omega"] / (1 - p["beta[1]"]))  # mean level
        e = eps / np.sqrt(s2)
        ln = (p["omega"]
              + p["alpha[1]"] * (np.abs(e) - np.sqrt(2 / np.pi))
              + p["gamma[1]"] * e + p["beta[1]"] * np.log(s2))
        return np.sqrt(252 * np.exp(ln))
    P = persistencia(f)
    s2 = p["omega"] / (1 - P)
    a = p["alpha[1]"] + p.get("gamma[1]", 0) * (eps < 0)
    return np.sqrt(252 * (p["omega"] + a * eps**2 + p["beta[1]"] * s2))

eps = np.array([-3.0, 0.0, 3.0])
pd.DataFrame({k: impacto(ajustes[k], eps)
              for k in ["GARCH-t", "GJR-t", "EGARCH-t"]},
             index=["3% drop", "no news", "3% rise"]).round(1)

         GARCH-t  GJR-t  EGARCH-t
3% drop     22.8   24.6      25.1
no news     17.5   17.4      18.9
3% rise     22.8   19.5      21.4

In the GARCH-t, a 3% drop and a 3% rise take volatility to the same place (22.8%): the curve is a symmetric parabola. In the GJR, the drop takes it to 24.6% and the rise to 19.5%, close to a day without news (17.4%): with $\alpha = 0.033$, rises barely move the variance. The EGARCH is asymmetric in the same way (25.1% and 21.4%), with a curve that grows exponentially rather than quadratically at the extremes. Lesson 1.6.2 found the same asymmetry without a model (tomorrow's standardized square was 1.84 after a drop and 1.28 after a rise); here it is quantified in a parameter.

Forecasting volatility out of sample

A volatility model is used to forecast, so it is evaluated by forecasting. The models are estimated with 2006-2018 and the variance of each session of 2019-2025 is forecast with fixed parameters (fit(last_obs=…) and forecast(start=…), which on each date uses only the information up to the previous day). The benchmarks are the EWMA of lesson 1.6.1 ($\lambda = 0.94$) and the variance of the previous 63 sessions.

The actual data point is $\varepsilon_t^2$, a very noisy proxy for that day's variance. With a noisy proxy, not every loss ranks the models the same way the true variance would; Patton (2011) derives which ones do, and QLIKE is one of them:

with $h_t$ the forecast variance (up to a constant that does not depend on the model). It depends on the ratio $\varepsilon_t^2/h_t$, so crisis days, with variances tens of times larger, do not dominate the mean as they do with the squared error. And it penalizes falling short more than overshooting: forecasting little volatility the day before a large move is expensive, which is what matters for risk.

import statsmodels.api as sm

corte = "2019-01-01"
prev = {}
for k, v in especs.items():
    f = arch_model(r, mean="Constant", **v).fit(last_obs=corte,
                                                disp="off")
    prev[k] = f.forecast(start=corte, reindex=False).variance["h.1"]
prev = pd.DataFrame(prev)
mu = r[:"2018"].mean()
e2 = (r - mu) ** 2
prev["EWMA"] = e2.ewm(alpha=0.06, adjust=False).mean().shift(1)
prev["63 sessions"] = e2.rolling(63).mean().shift(1)
real = e2[corte:]
prev = prev.reindex(real.index)

qlike = prev.apply(lambda h: real / h + np.log(h))
ecm = prev.apply(lambda h: (real - h) ** 2)
print(pd.DataFrame({"QLIKE": qlike.mean(), "MSE": ecm.mean()}).round(3))

def dm(a, b):  # Diebold-Mariano with Newey-West (lesson 2.1.3)
    d = (qlike[a] - qlike[b]).to_numpy()
    f = sm.OLS(d, np.ones(len(d))).fit(cov_type="HAC",
                                       cov_kwds={"maxlags": 10})
    return f"t = {f.tvalues[0]:.2f}"

print("GARCH-t versus EWMA:", dm("GARCH-t", "EWMA"))
print("GJR-t versus GARCH-t:", dm("GJR-t", "GARCH-t"))

             QLIKE     MSE
GARCH-N      0.849  16.496
GARCH-t      0.859  16.943
GJR-t        0.914  17.332
EGARCH-t     0.928  17.349
EWMA         1.100  22.017
63 sessions  1.150  23.781
GARCH-t versus EWMA: t = -6.21
GJR-t versus GARCH-t: t = 4.34

The four GARCH models clearly beat the EWMA and the 63-session window: the QLIKE loss drops from 1.10 to 0.86 with the GARCH-t, and the Diebold-Mariano test gives a t of −6.21. The difference lies in the long-run level: after a crisis, the EWMA is slow to come down and before another one it has nothing to return to.

The best in QLIKE is, in fact, the GARCH with normal innovations (0.849): to forecast tomorrow's variance, the distribution of the shocks matters little; it matters for the likelihood, the VaR and the tails.

The uncomfortable result is the other one: out of sample, the GJR forecasts worse than the symmetric GARCH (t = 4.3), even though in sample it improved on it by 26 points. Where does it lose?

dif = (qlike["GJR-t"] - qlike["GARCH-t"]).sort_values()
peores = dif.index[-4:]
pd.DataFrame({
    "rent_%": r[peores],
    "vol_GARCH_%": np.sqrt(252 * prev.loc[peores, "GARCH-t"]),
    "vol_GJR_%": np.sqrt(252 * prev.loc[peores, "GJR-t"]),
    "dif_QLIKE": dif[peores]}).round(2)

            rent_%  vol_GARCH_%  vol_GJR_%  dif_QLIKE
fecha
2025-02-04    2.95        23.75      14.92       4.95
2021-10-11    3.50        22.41      16.09       5.04
2019-09-04    3.24        20.47      12.57       9.32
2025-01-24    4.41        25.59      15.31      12.27

They are four large rises (from 2.95% to 4.41% in a day) for which the GJR forecast an annual volatility of 13% to 16% and the GARCH, of 20% to 26%: after bull streaks, the GJR, which hardly raises the variance on rises, expects very little risk. QLIKE penalizes that error heavily. On the course index, the asymmetry is real (it is planted and the test detects it), but it does not translate into better one-day forecasts; on the level 2 index futures (IDA, IDB, IDC), GARCH and GJR tie out of sample. Winning in sample is not enough: the comparison that counts is the out-of-sample one, with a suitable loss.

How long a volatility crisis lasts

The multi-day forecast comes from the persistence formula, and forecast(horizon=…) computes it. From the day with the highest conditional volatility in 2020:

gjr = ajustes["GJR-t"]
cv = gjr.conditional_volatility
pico = cv["2020"].idxmax()
f = gjr.forecast(start=pico, horizon=120, reindex=False)
v = f.variance.iloc[0]
p = gjr.params
lp = p["omega"] / (1 - persistencia(gjr))
mitad = int(np.argmax((v - lp) / (v.iloc[0] - lp) <= 0.5)) + 1
print(f"peak: {pico:%Y-%m-%d}, 1-day forecast vol "
      f"{np.sqrt(252 * v.iloc[0]):.1f}%; long run "
      f"{np.sqrt(252 * lp):.1f}%")
print(f"the deviation halves in {mitad} sessions")
np.sqrt(252 * v.iloc[[0, 9, 27, 59, 119]]).round(1)

peak: 2020-03-26, 1-day forecast vol 86.8%; long run 18.2%
the deviation halves in 30 sessions
h.001    86.8
h.010    78.1
h.028    63.5
h.060    44.8
h.120    26.7
Name: 2020-03-26 00:00:00, dtype: float64

On March 26, 2020, the GJR forecasts an annual volatility of 87% for the next day, and the forecast falls toward the long-run 18% with a half-life of 30 sessions: 63% at 28 sessions ahead, 27% at 120. For a system, that curve says how long a volatility-based sizing will keep the position reduced, and how much risk is expected over the horizon of the trade, not just tomorrow.

Volatility in a system

flowchart LR
  R[Returns] --> M[GARCH t]
  M --> S[Tomorrow's vol.]
  S --> Z[Size = target/vol]

Two frequent mistakes. The first, estimating the GARCH with the whole sample and using its conditional volatility in the sizing backtest: the parameters know the future (lesson 2.5.1); they have to be estimated with earlier data, as in the forecasting section. The second, applying it to a stock with earnings jumps: the jump on the announcement day is predictable in its date and the GARCH treats it as a shock that raises volatility for weeks (exercise 3).

Exercises

1. The fictional ETF KAPA (indice2_5min.parquet) has a daily GJR-GARCH-t volatility with $\alpha = 0.03$, $\gamma = 0.10$, $\beta = 0.86$ and $\nu = 6$ (level 2 README and generator), but only four years of data. Build its close-to-close daily returns, fit a GJR-t and compare the parameters with the truth using their standard errors.

k = pd.read_parquet("indice2_5min.parquet")
sesion = k["ts"].dt.tz_convert("America/New_York").dt.date
cierre = k.groupby(sesion)["cierre"].last()
rk = 100 * np.log(cierre).diff().dropna()
fk = arch_model(rk, mean="Constant", **especs["GJR-t"]).fit(
    disp="off")
pd.DataFrame({"estimado": fk.params, "error": fk.std_err,
              "verdad": [np.nan, np.nan, 0.03, 0.10, 0.86, 6]}
             ).round(3)

          estimado  error  verdad
mu           0.078  0.032     NaN
omega        0.115  0.054     NaN
alpha[1]     0.052  0.034    0.03
gamma[1]     0.030  0.037    0.10
beta[1]      0.832  0.063    0.86
nu           9.130  2.895    6.00

2. Scale the index to a target volatility of 10% a year with the conditional volatility of the GJR and with the EWMA, $w_t = 10 / \hat\sigma_t$, and measure the realized volatility of the scaled series overall and year by year. Which one comes closer to the target?

s2 = {"GJR": ajustes["GJR-t"].conditional_volatility ** 2,
      "EWMA": ((r - r.mean()) ** 2).ewm(alpha=0.06,
                                        adjust=False).mean().shift(1)}
for nombre, v in s2.items():
    esc = (10 / np.sqrt(252 * v) * r).dropna()
    anual = esc.groupby(esc.index.year).std() * np.sqrt(252)
    print(f"{nombre:5}: total {esc.std() * np.sqrt(252):.2f}%; "
          f"years from {anual.min():.1f} to {anual.max():.1f}%")
bruta = r.groupby(r.index.year).std() * np.sqrt(252)
print(f"unscaled: years from {bruta.min():.1f} to {bruta.max():.1f}%")

GJR  : total 10.04%; years from 8.6 to 11.5%
EWMA : total 10.69%; years from 9.8 to 12.8%
unscaled: years from 11.0 to 41.0%

3. In universo2 each stock has a jump on the day of the reaction to its earnings (level 2 README). What does a GARCH estimated on the returns of a stock do on the day after that jump, and how would you correct it?