---
title: "2.1.4 Volatility: ARCH, GARCH, EGARCH and GJR"
description: "GARCH(1,1) with t, GJR and EGARCH with the arch library: likelihood, persistence, half-life, news impact curve and out-of-sample QLIKE."
url: "https://pythonparatrading.com/en/python-trading-course/level-2/volatility-models-arch-garch-egarch-gjr/"
lang: "en"
image: "https://pythonparatrading.com/og/en/python-trading-course/level-2/volatility-models-arch-garch-egarch-gjr.png"
---

# 2.1.4 · Volatility models: ARCH, GARCH, EGARCH and GJR with the arch library

> Volatility changes, comes in streaks and rises more after drops. GARCH models turn those regularities into a forecast of tomorrow's variance. This lesson fits four variants to a synthetic index whose true volatility is known, measures how long shocks last and checks out of sample which model forecasts best.

- Course: [Level 2: Systematic strategies and rigorous backtesting](https://pythonparatrading.com/en/python-trading-course/level-2/)
- Module: 2.1 · Statistics and econometrics for trading
- Duration: 70 min
- Prerequisites: [1.6.1 · Returns, volatility, rolling correlations and drawdown](https://pythonparatrading.com/en/python-trading-course/level-1/returns-volatility-correlation-and-drawdown/), [1.6.2 · Stylized facts: fat tails, clustering, autocorrelation and leverage](https://pythonparatrading.com/en/python-trading-course/level-1/stylized-facts-of-financial-markets/), [2.1.1 · Hypothesis tests, p-values, bootstrap and regression with robust standard errors](https://pythonparatrading.com/en/python-trading-course/level-2/hypothesis-testing-bootstrap-robust-errors/)

## By the end of this lesson

- Write and estimate with `arch` a GARCH(1,1) with normal and Student's t innovations, a GJR-GARCH and an EGARCH, and compare them by likelihood, AIC and BIC.
- Compute the persistence, the long-run variance and the half-life of a volatility shock, and explain why regime changes inflate it.
- Draw the news impact curve of each model and relate it to the leverage effect.
- Evaluate out-of-sample volatility forecasts with the QLIKE loss against the EWMA and a rolling window.

## 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$):

$$
\sigma_t^2 = \omega + \alpha\,\varepsilon_{t-1}^2 + \beta\,\sigma_{t-1}^2
$$

- $\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:

```mermaid
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)[^1]: 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$:

$$
\ln\sigma_t^2 = \omega + \alpha\left(|e_{t-1}| - \sqrt{2/\pi}\right) + \gamma\, e_{t-1} + \beta\ln\sigma_{t-1}^2
$$

  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.

```python
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}")
```
```output
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).

```python
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)
```
```output
          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.

```python
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}")
```
```output
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

$$
\mathrm{E}_t\left[\sigma_{t+h}^2\right] = \sigma^2_{LP} + P^{\,h}\left(\sigma_t^2 - \sigma^2_{LP}\right)
$$

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$.

```python
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)
```
```output
                    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):

```python
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)
```
```output
                         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:

*Fig. — Annual conditional volatility of the GJR-t on the synthetic equal-weighted index (every 10 sessions) against the market factor's long-run volatility in each segment of the generator's script, 2006-2025.* A jagged conditional volatility line and a step line with the script's long-run volatility. The step goes from 12% in calm periods to 52% in the 2008 crisis, 36% in 2011 and 80% in March 2020, the highlighted segment. The conditional volatility follows the steps with a lag and its own spikes: it reaches 90% in March 2020 and comes down slowly afterwards.

## 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.

```python
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)
```
```output
         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:

$$
L(\varepsilon_t^2, h_t) = \frac{\varepsilon_t^2}{h_t} + \ln h_t
$$

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.

```python
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"))
```
```output
             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
```

*Fig. — Mean out-of-sample QLIKE loss (2019-2025) of six forecasts of the synthetic equal-weighted index's daily variance; lower is better.* Horizontal bars: GARCH with normal 0.849 (the best, highlighted), GARCH-t 0.859, GJR-t 0.914, EGARCH-t 0.928, EWMA 1.100 and 63-session variance 1.150. The four GARCH models sit well below the EWMA and the rolling window.

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?

```python
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)
```
```output
            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:

```python
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)
```
```output
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
```

*Fig. — Annual volatility forecast by the GJR-t from the 26 March 2020 peak for each of the next 120 days, against its long-run level (synthetic equal-weighted index).* Declining curve with the horizon from 1 to 120 days: it starts at 86.8% annual volatility, passes 62% at 30 days (highlighted: the variance's deviation from the long run has halved), 44.8% at 60 and 26.7% at 120, approaching a dashed line at 18.2%, the long-run volatility.

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

| Use | What a GARCH contributes | Caveat |
|---|---|---|
| Volatility-based *sizing* (lesson 1.6.1) | Forecast with mean reversion: it does not stay low after a crisis | Re-estimate with a long window; the parameters change little |
| VaR and stops | Tomorrow's volatility and a t distribution with its ν | Tails: check against the worst real days (lesson 2.5.7) |
| Multi-day horizon | Term structure of the forecast | Persistence includes the regimes: inflated half-life |
| Synthetic data for stress tests | Simulate series with clustering and asymmetry (lesson 2.5.7) | The parameters of one stretch are not valid for another |
| Options | Compare forecast volatility with implied volatility (level 3) | GARCH is a physical forecast, not a risk-neutral one |

```mermaid
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.

**Solution.**

```python
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)
```
```output
          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
```

With 1,002 sessions, the standard errors are the size of the parameters themselves: $\gamma$ comes out at 0.030 with
an error of 0.037, so the 95% interval runs from −0.04 to 0.10 and rules out neither symmetry nor the truth. $\nu$
comes out at 9.1 with an error of 2.9. A GARCH needs many years to estimate the asymmetry and the tails well; with
short series, it is better to fix reasonable parameters or to estimate them with a similar asset with a long history.

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?

**Solution.**

```python
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}%")
```
```output
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%
```

Both stabilize the risk a lot: unscaled, the annual volatility of the index ranges from 11% to 41%; scaled, from 8.6%
to 11.5% with the GJR and from 9.8% to 12.8% with the EWMA. The GJR hits the level better (10.04% versus 10.69%),
because the EWMA falls short at the onset of crises. Note: here the GJR is estimated with the whole sample, so its
advantage is somewhat inflated; in a real backtest it has to be estimated with the past only.

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?

**Solution.**

The jump enters as an enormous $\varepsilon^2$ and the GARCH, which does not know it was an earnings day, raises the
forecast variance for the following days and lets it come down with its half-life: it overestimates the risk for
weeks and a volatility-based sizing reduces the position for no reason. In addition, the jumps inflate $\alpha$ and
the tails of the t. The usual corrections are three: add to the variance equation a variable equal to 1 on the
earnings day (`arch_model` accepts regressors `x` in the mean, not in the variance, so that extension has to be coded),
exclude those days from the estimation and treat their risk separately (lesson 2.3.5), or estimate the GARCH on the
index or the sector, without idiosyncratic jumps, and add the event risk only when it is due.

[^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>.

## Self-assessment

1. A GARCH(1,1) has α = 0.10 and β = 0.88. What is the half-life of a volatility shock?

   - A. About 2 sessions
   - B. About 98 sessions
   - C. About 34 sessions: the deviation of the variance from its long-run level is multiplied by α + β = 0.98 every day, and ln 0.5 / ln 0.98 ≈ 34
   - D. Infinite

   *Answer C: About 34 sessions: the deviation of the variance from its long-run level is multiplied by α + β = 0.98 every day, and ln 0.5 / ln 0.98 ≈ 34.* The h-day forecast of the variance is the long-run one plus (α + β) to the power of h times today's deviation from it. The half-life is the h that leaves the deviation at half. With persistence 0.98, 34 sessions; with the true persistence of the market factor in the data (0.945), 12.

2. The GJR estimated on the equal-weighted index has a persistence of 0.976, but the generator's market factor has 0.945. What explains the difference?

   - A. A bug in the arch library
   - B. That the generator's long-run volatility changes in stretches (calm, crisis): a GARCH with a single long-run level interprets those regime changes as very persistent shocks
   - C. That the index has more stocks
   - D. That Student's t has fat tails

   *Answer B: That the generator's long-run volatility changes in stretches (calm, crisis): a GARCH with a single long-run level interprets those regime changes as very persistent shocks.* It is a known result (Lamoureux and Lastrapes, 1990). In a calm stretch with constant long-run volatility, the estimated persistence drops to 0.96. A half-life of 28 sessions mixes the dynamics of the shocks with that of the regimes.

3. In the GJR news impact curve, a 3% drop takes tomorrow's annual volatility to 24.6% and a 3% rise, to 19.5%. Which parameter produces the difference?

   - A. β
   - B. ω
   - C. The degrees of freedom of the t
   - D. γ: the term that only acts when the shock is negative

   *Answer D: γ: the term that only acts when the shock is negative.* In the GJR, the effect of ε² on tomorrow's variance is α if ε > 0 and α + γ if ε < 0. With α = 0.033 and γ = 0.100, a drop weighs four times as much as a rise of the same size: it is the leverage effect that lesson 1.6.2 measured without a model.

4. Why is the QLIKE loss used to compare variance forecasts instead of the squared error?

   - A. Because the actual data point (r²) is a very noisy proxy for the variance: QLIKE ranks the models the same as if the true variance were known and, since it depends on the ratio r²/h, it does not let a few crisis days dominate the mean
   - B. Because it is easier to compute
   - C. Because it always favors the GARCH
   - D. Because it does not need the return

   *Answer A: Because the actual data point (r²) is a very noisy proxy for the variance: QLIKE ranks the models the same as if the true variance were known and, since it depends on the ratio r²/h, it does not let a few crisis days dominate the mean.* Patton (2011) derives which losses rank models the same as with the true variance even though the proxy is noisy; QLIKE is one of them. Since it depends on r²/h, a day with a variance a hundred times larger does not weigh a hundred times more, and it penalizes falling short more than overshooting.

5. In sample, the GJR improves the GARCH log-likelihood by 26 points; out of sample (2019-2025), its QLIKE is worse. Which conclusion is correct?

   - A. That the GJR is badly implemented
   - B. That the leverage effect does not exist in these data
   - C. That the EWMA should be used
   - D. That the best in-sample fit does not guarantee a better forecast: on this index, the GJR forecasts very little volatility after bull streaks and pays for it on days with large unexpected moves. The honest comparison is out of sample

   *Answer D: That the best in-sample fit does not guarantee a better forecast: on this index, the GJR forecasts very little volatility after bull streaks and pays for it on days with large unexpected moves. The honest comparison is out of sample.* The asymmetry is real (it is planted and the likelihood ratio test detects it), but forecasting is a different task. On the course's index futures, GARCH and GJR tie out of sample; both clearly beat the EWMA and the 63-session window.

## Key points

- GARCH(1,1): today's variance = ω + α · (yesterday's shock)² + β · (yesterday's variance). α measures the reaction to news; β, the memory; α + β, the persistence.
- The EWMA of lesson 1.6.1 is a GARCH without ω and without a long-run level: it does not return to any mean. The GARCH does, to ω/(1 − α − β).
- Student's t improves the log-likelihood by 185 points versus the normal: fat tails do not disappear when standardizing.
- GJR and EGARCH capture the leverage effect; on the course index, γ = 0.100 versus the generator's 0.11.
- The estimated persistence (0.976, half-life of 28 sessions) is inflated by the regime changes in long-run volatility; in a calm stretch it drops to 0.96.
- Out of sample, the GARCH models reduce the QLIKE loss versus the EWMA and the rolling window; asymmetry does not always help forecasting.

## Resources

### Research

- [Volatility Estimators for Trading](https://pythonparatrading.com/en/volatility-estimators-for-trading/) · Research-0002 · 4 min

### Documentation

- [arch: univariate volatility models (documentation)](https://arch.readthedocs.io/en/latest/univariate/introduction.html) · arch.readthedocs.io

### Books

- [Bollerslev, T. (1986): Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 31(3), 307-327](https://doi.org/10.1016/0304-4076(86)90063-1) · doi.org
- [Glosten, L. R., Jagannathan, R. and Runkle, D. E. (1993): On the Relation between the Expected Value and the Volatility of the Nominal Excess Return on Stocks. Journal of Finance 48(5), 1779-1801](https://doi.org/10.1111/j.1540-6261.1993.tb05128.x) · doi.org
- [Lamoureux, C. G. and Lastrapes, W. D. (1990): Persistence in Variance, Structural Change, and the GARCH Model. Journal of Business & Economic Statistics 8(2), 225-234](https://doi.org/10.1080/07350015.1990.10509794) · doi.org
- [Patton, A. J. (2011): Volatility forecast comparison using imperfect volatility proxies. Journal of Econometrics 160(1), 246-256](https://doi.org/10.1016/j.jeconom.2010.03.034) · doi.org

### Data

- [Level 2 data: adjusted prices of the synthetic universe of 120 stocks (Parquet)](https://pythonparatrading.com/curso/datos/nivel-2/universo2_ajustados.parquet) · Download · universo2_ajustados.parquet

Previous lesson: [2.1.3 · ARMA and ARIMA models: autocorrelation (ACF) and partial autocorrelation (PACF)](https://pythonparatrading.com/en/python-trading-course/level-2/arma-arima-models-acf-pacf/)

Next lesson: [2.1.5 · Cointegration: Engle-Granger, Johansen and the half-life of mean reversion](https://pythonparatrading.com/en/python-trading-course/level-2/cointegration-engle-granger-johansen-half-life/)