Asymmetrische en zwaarstaartige GARCH: EGARCH, GJR en Student-t

In Deel 1 van deze reeks hebben we GARCH(1,1) van de grond af opgebouwd: de intuïtie achter volatiliteitsclustering, de recursie van de conditionele variantie, maximum likelihood-schatting, forecasting en de standaard residudiagnostiek met de arch-bibliotheek. Als je dat nog niet gelezen hebt, begin daar — deze post gaat ervan uit dat je een gewone GARCH(1,1) al kunt fitten en interpreteren en herleidt de basis niet opnieuw.

Gewone GARCH(1,1) is een goede baseline en een slecht eindantwoord. Het heeft twee structurele gebreken die goedkoop te negeren zijn in een backtest en duur te negeren met echt kapitaal. Ten eerste is het symmetrisch: het model reageert op een $+5\%$-dag precies zoals het reageert op een $-5\%$-dag, omdat de schok alleen via zijn kwadraat $\varepsilon_{t-1}^2$ in de variantierecursie terechtkomt. Kwadrateren gooit het teken weg. Ten tweede gaat het uit van Gaussische innovaties: zelfs nadat GARCH de volatiliteitsclustering heeft geabsorbeerd, zijn de gestandaardiseerde residuen van BTC en ETH zichtbaar dikstaartig, en een Gaussische likelihood onderschat de staart systematisch. Een GARCH(1,1)-Normal 99%-VaR wordt veel vaker doorbroken dan 1% van de tijd.

Deze post lost beide gebreken op. We voegen asymmetrie toe met GJR-GARCH en EGARCH, en dikke staarten met Student-$t$ en Hansens scheve-$t$ innovaties. Vervolgens doen we het ding dat daadwerkelijk de huur betaalt: de gefitte conditionele verdeling omzetten in een eenstaps-Value-at-Risk- en Expected Shortfall-forecast, en die forecast eerlijk backtesten met de Kupiec- en Christoffersen-toetsen. Een volatiliteitsmodel dat je nooit risicotest is een versiersel.

Het hefboomeffect, en waarom crypto rommeliger is

In aandelen heeft de asymmetrie een naam en een verhaal. Het hefboomeffect (Black, 1976): wanneer het aandeel van een bedrijf daalt, stijgt zijn schuld/eigen-vermogen-verhouding, wordt het eigen vermogen mechanisch risicovoller en neemt de volatiliteit toe. Slecht nieuws verhoogt de toekomstige volatiliteit meer dan even groot goed nieuws. Empirisch is dit een van de robuustste gestileerde feiten in de literatuur over aandelenvolatiliteit.

Crypto heeft geen eigen vermogen en geen balanshefboom in de bedrijfsmatige zin, en toch vertoont zich meestal nog steeds een hefboomeffect-achtige asymmetrie — gedreven door gedwongen deleveraging in plaats van boekhouding. Wanneer BTC hard daalt, worden overgecollateraliseerde leningen geliquideerd, worden long-posities in perpetual futures gedwongen gesloten, klapt de funding om, en voedt de cascade de volatiliteit. Het mechanisme verschilt dus, maar het teken stemt vaak overeen met aandelen: neerwaartse bewegingen jagen de volatiliteit sterker op.

De belangrijke kanttekening: crypto is rommeliger, en je zou asymmetrie moeten behandelen als een empirische vraag in plaats van een wet. Heftige opwaartse bewegingen — short squeezes, een door hefboom aangedreven melt-up, een gap door de goedkeuring van een ETF — kunnen ook de gerealiseerde volatiliteit opjagen. Afhankelijk van het activum en het steekproefvenster kan de geschatte asymmetrie sterk, zwak of af en toe het "verkeerde" teken hebben. De discipline waar deze post op staat: fit het asymmetrische model, kijk of de asymmetrieparameter statistisch significant is en in de verwachte richting wijst, en houd de extra parameter alleen aan als hij zijn plaats verdient. Ga er niet van uit dat het aandelenverhaal overdraagbaar is; toets het.

Toetsen op asymmetrie voordat je het modelleert

De briefing hierboven zegt "behandel asymmetrie als empirisch" — dus voordat je een asymmetrisch model fit, voer je een goedkope formele toets uit op of asymmetrie überhaupt aanwezig is. De Engle-Ng sign-bias-toetsen (1993) doen precies dit. Fit eerst een symmetrische GARCH(1,1), neem de gekwadrateerde gestandaardiseerde residuen $z_t^2$, en regresseer ze op indicatoren van het teken en de grootte van de vorige schok:

$z_t^2 = a_0 + a_1 S_{t-1}^- + a_2 S_{t-1}^- \varepsilon_{t-1} + a_3 S_{t-1}^+ \varepsilon_{t-1} + u_t$

waarbij $S_{t-1}^- = \mathbf{1}\{\varepsilon_{t-1} < 0\}$ en $S_{t-1}^+ = 1 - S_{t-1}^-$. De logica: als het symmetrische model alles al heeft opgevangen, mogen het teken en de grootte van de schok van gisteren het gekwadrateerde residu van vandaag niet voorspellen, dus $a_1 = a_2 = a_3 = 0$. Individuele $t$-toetsen zijn de sign-bias ($a_1$), negative-size-bias ($a_2$) en positive-size-bias ($a_3$) toetsen; een gezamenlijke $F$-toets op alle drie is de omnibus. Een significante $a_1$ of $a_2$ zegt dat negatieve schokken systematisch verkeerd worden geprijsd door het symmetrische model — je aanwijzing dat GJR of EGARCH zal helpen.

import statsmodels.api as sm

def sign_bias_test(symmetric_res):
    """Engle-Ng sign-bias tests on a fitted symmetric GARCH result."""
    z = symmetric_res.std_resid.dropna()
    z2 = (z ** 2).values[1:]
    eps_lag = z.values[:-1]                      # standardized shock proxy
    neg = (eps_lag < 0).astype(float)
    X = np.column_stack([
        np.ones_like(eps_lag),                   # intercept
        neg,                                     # sign bias
        neg * eps_lag,                           # negative size bias
        (1 - neg) * eps_lag,                     # positive size bias
    ])
    ols = sm.OLS(z2, X).fit()
    names = ["const", "sign_bias", "neg_size_bias", "pos_size_bias"]
    for nm, coef, t, p in zip(names, ols.params, ols.tvalues, ols.pvalues):
        print(f"{nm:16s} coef={coef:+.4f}  t={t:+.2f}  p={p:.3f}")
    print(f"Joint F p-value: {ols.f_pvalue:.4f}")
    return ols

sign_bias_test(models["GARCH-N"])

Als de gezamenlijke $F$-toets niet significant is, heb je empirische toestemming om symmetrisch te blijven en twee parameters te besparen. Als hij wel significant is — de gebruikelijke uitkomst voor BTC/ETH — ga dan met een gerust geweten door naar GJR/EGARCH, wetende dat je een echt kenmerk modelleert en geen ruis najaagt. Dit is de empirische discipline die de aanhef eiste: ga niet uit van het hefboomverhaal uit de aandelenwereld, maar toets ernaar.

GJR-GARCH: asymmetrie via een drempelterm

Het Glosten-Jagannathan-Runkle-model (1993) — soms TGARCH of threshold GARCH genoemd — is de kleinst mogelijke aanpassing aan GARCH(1,1) die slecht nieuws en goed nieuws verschillende effecten laat hebben. Herinner je de symmetrische recursie van de conditionele variantie uit Deel 1:

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

GJR voegt één drempelterm toe: een extra dosis variantie die alleen inschakelt na een negatieve schok.

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

waarbij $I_{t-1}$ de indicator is

$I_{t-1} = \begin{cases} 1 & \text{if } \varepsilon_{t-1} < 0 \\ 0 & \text{if } \varepsilon_{t-1} \geq 0 \end{cases}$

Lees de recursie per geval. Na een positieve schok ($\varepsilon_{t-1} \geq 0$) is de indicator nul en is de impact van de gekwadrateerde schok op de variantie van de volgende periode gewoon $\alpha$. Na een negatieve schok is de indicator één en is de impact $\alpha + \gamma$. De parameter $\gamma$ is het hele asymmetrieverhaal in één getal:

• $\gamma > 0$: negatieve schokken verhogen de volatiliteit meer dan positieve schokken van dezelfde grootte. Dit is het hefboomeffect, en het is wat je meestal verwacht te vinden in BTC/ETH.
• $\gamma = 0$: het model valt terug op symmetrische GARCH(1,1). Een likelihood-ratio- of $t$-toets op $\gamma$ is daarom een directe toets op of asymmetrie überhaupt bestaat.
• $\gamma < 0$: positieve schokken verhogen de volatiliteit meer — het incidentele crypto-melt-up-regime. Zeldzaam, maar sluit het niet a priori uit.

Positiviteit en stationariteit

Omdat $\sigma_t^2$ nog steeds additief wordt opgebouwd, hebben we nodig dat elke term niet-negatief blijft. De voldoende positiviteitsvoorwaarden zijn

$\omega > 0, \quad \alpha \geq 0, \quad \alpha + \gamma \geq 0, \quad \beta \geq 0.$

Merk op dat $\gamma$ zelf negatief mag zijn zolang $\alpha + \gamma \geq 0$, zodat de impact na slecht nieuws nooit negatief wordt.

Voor covariantie-stationariteit nemen we aan dat de innovaties $z_t = \varepsilon_t/\sigma_t$ gestandaardiseerd zijn met een symmetrische verdeling rond nul, zodat $P(\varepsilon_{t-1} < 0) = 1/2$ en de indicator gemiddeld $\gamma/2$ bijdraagt. De stationariteitsvoorwaarde wordt

$\alpha + \beta + \tfrac{1}{2}\gamma < 1.$

De onconditionele (langetermijn-)variantie is dan

$\bar{\sigma}^2 = \frac{\omega}{1 - \alpha - \beta - \tfrac{1}{2}\gamma}.$

Dit is de GJR-analogon van het resultaat $\bar{\sigma}^2 = \omega/(1-\alpha-\beta)$ uit Deel 1, met de extra $\tfrac{1}{2}\gamma$-term die de gemiddelde bijdrage van de hefboom-halfwaardetijd verantwoordt. Als je innovatieverdeling scheef is (Hansens skew-$t$, hieronder), wordt de $1/2$ vervangen door de daadwerkelijke kans dat $z_t < 0$, maar $1/2$ is de standaardreferentie die wordt gebruikt voor de gerapporteerde persistentie.

EGARCH: log-variantie modelleren, geen positiviteitsbeperkingen

GJR houdt je in een dwangbuis van variantie-positiviteit: elke parametercombinatie moet worden gecontroleerd tegen ongelijkheidsbeperkingen, wat lastig is tijdens optimalisatie en erger tijdens rollende herschatting wanneer een venster af en toe in een onhaalbaar gebied afdwaalt. Nelsons Exponential GARCH (1991) omzeilt dit volledig door de logaritme van de conditionele variantie te modelleren. Omdat $\log \sigma_t^2$ elk reëel getal kan zijn, is $\sigma_t^2 = \exp(\log \sigma_t^2)$ automatisch positief ongeacht de parameters. Geen beperkingen om op te leggen.

Schrijf de recursie in termen van de gestandaardiseerde innovatie $z_{t-1} = \varepsilon_{t-1}/\sigma_{t-1}$:

$\log \sigma_t^2 = \omega + \beta \log \sigma_{t-1}^2 + \alpha\bigl(|z_{t-1}| - \mathbb{E}|z_{t-1}|\bigr) + \gamma\, z_{t-1}$

Twee termen dragen de schok, en het scheiden ervan is het hele idee:

• De magnitude-term $\alpha(|z_{t-1}| - \mathbb{E}|z_{t-1}|)$ reageert op de grootte van de schok, met het teken verwijderd. Aftrekken van $\mathbb{E}|z_{t-1}|$ centreert hem zodat een schok van gemiddelde grootte niets bijdraagt. Voor een standaardnormale verdeling is $\mathbb{E}|z| = \sqrt{2/\pi} \approx 0.7979$; voor een gestandaardiseerde Student-$t$ is de verwachte absolute waarde kleiner en hangt af van $\nu$, maar arch handelt dit intern af.
• De teken-term $\gamma\, z_{t-1}$ is de asymmetrie. Hij is lineair in de innovatie met teken, dus een negatieve $z_{t-1}$ duwt $\log\sigma_t^2$ in de tegengestelde richting van een positieve.

De tekenconventie is van belang en brengt mensen in de war. In deze parametrisatie komt het hefboomeffect (slecht nieuws verhoogt de volatiliteit) overeen met $\gamma < 0$: een negatieve schok $z_{t-1} < 0$ maakt dan $\gamma z_{t-1} > 0$, wat de log-variantie verhoogt. Dit is het tegenovergestelde teken van GJR's $\gamma > 0$. Lees altijd de eigen documentatie van het model voor de conventie in plaats van aan te nemen; arch rapporteert EGARCH met zijn eigen teken, en we controleren het hieronder tegen een news impact curve in plaats van op ons geheugen te vertrouwen.

Omdat alles additief is in logs, wordt de persistentie van een EGARCH(1,1) bepaald door de enkele autoregressieve coëfficiënt $\beta$ op $\log\sigma_{t-1}^2$; stationariteit vereist alleen $|\beta| < 1$. Dat is een veel schonere voorwaarde dan de GJR-ongelijkheid, en het is een echt praktisch voordeel wanneer je herfit op rollende vensters.

Een subtiliteit die het waard is om te noemen: de reactie van EGARCH op schokken is exponentieel in de innovatie (je exponentieert aan het einde), terwijl GJR kwadratisch is. EGARCH reageert daarom heftiger op grote schokken — een kenmerk in crypto, waar de staartgebeurtenissen de gebeurtenissen zijn die ertoe doen, maar ook een reden waarom EGARCH af en toe onwaarschijnlijk grote variantieforecasts kan produceren na een uitschieterdag. Geen van beide is universeel beter; je kiest op basis van out-of-sample fit en risicobacktests, wat het punt is van deze hele reeks.

De news impact curve

De schoonste manier om het verschil tussen symmetrische GARCH, GJR en EGARCH te zien is de news impact curve (Engle en Ng, 1993): houd $\sigma_{t-1}$ vast op zijn langetermijnniveau en plot de conditionele variantie van de volgende periode $\sigma_t^2$ als functie van de laatste schok $\varepsilon_{t-1}$. Hij beantwoordt "gegeven een schok van deze grootte en dit teken, hoeveel verhoogt het model de volatiliteit van morgen?"

• Symmetrische GARCH produceert een symmetrische parabool gecentreerd op nul. Een schok van $-5\%$ en een van $+5\%$ landen op dezelfde hoogte. Dit is precies het gebrek dat we oplossen.
• GJR produceert een parabool met een knik bij nul — steiler aan de linkerkant (negatieve schokken) dan aan de rechterkant wanneer $\gamma > 0$. De twee helften hebben respectievelijk kromming $\alpha+\gamma$ en $\alpha$.
• EGARCH produceert een asymmetrische, exponentiële V-vorm: de twee armen hebben verschillende hellingen vanwege de $\gamma z$-term, en het geheel kromt sneller omhoog dan een parabool vanwege de uiteindelijke exponentiatie.

We plotten alle drie op basis van gefitte parameters later, in de implementatiesectie — het is de nuttigste enkele diagnostiek om te communiceren wat asymmetrie je oplevert.

Dikke staarten: Student-t en scheve-t innovaties

Asymmetrie corrigeert de reactie van het model op het teken van schokken. Het doet niets aan de verdeling van de schokken zelf. Gewone GARCH gaat uit van $z_t \sim \mathcal{N}(0,1)$, en die aanname is bijna altijd verkeerd voor crypto. Zelfs nadat GARCH de volatiliteitsclustering heeft verwijderd, behouden de gestandaardiseerde residuen $z_t = \varepsilon_t / \sigma_t$ overtollige kurtosis — ze zijn dikstaartig. Een Gaussische likelihood, die de schouders van de verdeling fit, onderweegt hoe vaak een gestandaardiseerde dag van $4$, $5$ of $6$ sigma daadwerkelijk voorkomt.

Het gevolg voor risico is direct. Een Gaussische 99%-VaR gebruikt het kwantiel $\Phi^{-1}(0.01) \approx -2.326$, dus voorspelt hij $\text{VaR}_{0.99} \approx 2.326\,\sigma_t$. Als de ware gestandaardiseerde verdeling Student-$t$ is met bijvoorbeeld $\nu = 5$ vrijheidsgraden, ligt het ware 1%-kwantiel rond $-3.36$ — de Gaussische VaR is ongeveer $44\%$ te optimistisch bij dat betrouwbaarheidsniveau. Je zult hem veel vaker doorbreken dan 1% van de tijd en systematisch verrast worden door "onmogelijke" dagen. Dit is geen crypto-eigenaardigheid; Bollerslev (1987) introduceerde $t$-GARCH juist omdat aandelen- en FX-residuen dezelfde dikke staarten vertoonden. Crypto is gewoon een extremere versie van hetzelfde probleem.

Gestandaardiseerde Student-t

De Student-$t$-dichtheid heeft een vrijheidsgraden-parameter $\nu > 2$ die de staartdikte regelt: kleine $\nu$ betekent dikke staarten, en naarmate $\nu \to \infty$ convergeert de $t$ naar de Gaussische. De valkuil is dat de ruwe $t_\nu$-verdeling variantie $\nu/(\nu-2) \neq 1$ heeft, dus moeten we hem standaardiseren naar eenheidsvariantie voordat we hem als innovatie gebruiken — anders zou de "$\sigma_t$" in de GARCH-recursie niet daadwerkelijk de conditionele standaarddeviatie zijn.

De gestandaardiseerde Student-$t$-innovatie met eenheidsvariantie heeft dichtheid

$f(z;\nu) = \frac{\Gamma\!\left(\frac{\nu+1}{2}\right)}{\sqrt{\pi(\nu-2)}\;\Gamma\!\left(\frac{\nu}{2}\right)}\left(1 + \frac{z^2}{\nu-2}\right)^{-\frac{\nu+1}{2}}, \qquad \nu > 2.$

Let op de $(\nu-2)$ erin — dat is de standaardisatie, herschaald zodat $\mathrm{Var}(z)=1$. De bijdrage aan de log-likelihood van één observatie, gegeven de GARCH-conditionele variantie $\sigma_t^2$ en $z_t = \varepsilon_t/\sigma_t$, is

$\ell_t = \log\Gamma\!\left(\tfrac{\nu+1}{2}\right) - \log\Gamma\!\left(\tfrac{\nu}{2}\right) - \tfrac{1}{2}\log\bigl(\pi(\nu-2)\bigr) - \tfrac{1}{2}\log\sigma_t^2 - \frac{\nu+1}{2}\log\!\left(1 + \frac{z_t^2}{\nu-2}\right).$

De term $-\tfrac{1}{2}\log\sigma_t^2$ is de Jacobiaan van de transformatie van $\varepsilon_t$ naar $z_t$ — dezelfde term die je zag in de Gaussische GARCH-likelihood in Deel 1. Alleen de vorm verandert. Het maximaliseren van $\sum_t \ell_t$ gezamenlijk over de GARCH-parameters en $\nu$ is precies wat arch doet wanneer je dist='t' doorgeeft.

De geschatte $\nu$ is zelf informatief. Voor dagelijkse BTC/ETH-rendementen kom je doorgaans in het bereik $\nu \approx 3\text{–}6$ terecht — dikke staarten, maar met eindige variantie (waarvoor $\nu > 2$ nodig is) en meestal eindige kurtosis (waarvoor $\nu > 4$ nodig is). Als je gefitte $\nu$ onder de 4 zakt, wees je er dan van bewust dat de steekproefkurtosis in het model technisch gezien oneindig is en sommige schatters instabiel worden; het is een signaal om goed te kijken naar uitschieters en datakwaliteit.

Hansens scheve-t

Student-$t$ is dikstaartig maar nog steeds symmetrisch — de linker- en rechterstaart zijn even zwaar. Crypto-rendementsresiduen zijn vaak ook scheef: de linkerstaart (crashes) is zwaarder dan de rechter. Hansens scheve-$t$ (1994) veralgemeent de gestandaardiseerde $t$ met een scheefheidsparameter $\lambda \in (-1, 1)$ naast $\nu$:

$f(z;\nu,\lambda) = \begin{cases} b\,c\left(1 + \dfrac{1}{\nu-2}\left(\dfrac{bz+a}{1-\lambda}\right)^2\right)^{-\frac{\nu+1}{2}} & z < -a/b \\[2mm] b\,c\left(1 + \dfrac{1}{\nu-2}\left(\dfrac{bz+a}{1+\lambda}\right)^2\right)^{-\frac{\nu+1}{2}} & z \geq -a/b \end{cases}$

waarbij de constanten $a = 4\lambda c\,\frac{\nu-2}{\nu-1}$, $b^2 = 1 + 3\lambda^2 - a^2$ en $c = \frac{\Gamma(\frac{\nu+1}{2})}{\sqrt{\pi(\nu-2)}\,\Gamma(\frac{\nu}{2})}$ zo gekozen zijn dat $z$ nulgemiddelde en eenheidsvariantie heeft voor elke geldige $(\nu,\lambda)$. De verdeling splitst bij $z = -a/b$ en gebruikt in elk stuk een andere schaling om meer massa in de ene staart te buigen.

Interpretatie: $\lambda < 0$ geeft een linksscheve verdeling (zwaardere onderkant), wat de gebruikelijke bevinding voor crypto is en wat je zou verwachten te koppelen aan een hefboomeffect. $\lambda = 0$ herstelt de symmetrische Student-$t$, dus een toets op $\lambda = 0$ vertelt je of de scheefheidsterm iets oplevert. In arch is dit dist='skewt', dat zowel $\nu$ als $\lambda$ schat. De beloning is een VaR waarvan het linkerstaart-kwantiel eerlijk zwaarder is dan het rechterstaart-kwantiel — precies wat je wilt wanneer de verliezen die je probeert te overleven asymmetrisch zijn. Dit sluit direct aan bij de asymmetrie van verlies versus winst in positie-uitkomsten: een drawdown van $x\%$ heeft meer dan $x\%$ nodig om te herstellen, dus het verkeerd modelleren van de linkerstaart is kostbaarder dan het verkeerd modelleren van de rechter.

Python-implementatie

We fitten dit alles nu met de arch-bibliotheek. De opzet weerspiegelt Deel 1: haal dagelijkse rendementen op, schaal met 100 voor numerieke conditionering (GARCH-optimizers gedragen zich slecht wanneer rendementen $O(0.01)$ zijn), en fit met een constant gemiddelde. Als je intraday of een ander gemiddelde-model wilt, is de machinerie identiek.

Opzet en data

import numpy as np
import pandas as pd
from arch import arch_model
from arch.univariate import GARCH, EGARCH, ConstantMean, StudentsT, SkewStudent
from scipy import stats

def fetch_returns(symbol="BTC-USD", start="2019-01-01", end="2025-12-31"):
    """Daily log returns, in percent (scaled x100 for the optimizer)."""
    import yfinance as yf
    px = yf.download(symbol, start=start, end=end)["Close"].dropna()
    ret = 100.0 * np.log(px / px.shift(1)).dropna()
    ret.name = symbol
    return ret

r = fetch_returns("BTC-USD")
print(f"{len(r)} daily observations, "
      f"annualized vol ~ {r.std() * np.sqrt(365):.0f}%")

Crypto wordt 24/7 verhandeld, dus we annualiseren met 365, niet met 252 — een kleine maar terugkerende bron van verwarring wanneer je een crypto-Sharpe of -volatiliteit vergelijkt met de cijfers van een aandelen-desk.

Vier modellen fitten

Het patroon in arch: vol='Garch' met p=1, q=1 is symmetrische GARCH; toevoegen van o=1 schakelt de GJR-drempelterm in; vol='EGARCH' schakelt over naar het log-variantiemodel. De innovatieverdeling wordt ingesteld met dist: 'normal', 't', 'skewt'.

def fit(returns, vol="Garch", p=1, o=0, q=1, dist="normal"):
    am = arch_model(returns, mean="Constant",
                    vol=vol, p=p, o=o, q=q, dist=dist)
    res = am.fit(disp="off")
    return res

models = {
    "GARCH-N":    fit(r, vol="Garch",  o=0, dist="normal"),  # Part 1 baseline
    "GJR-t":      fit(r, vol="Garch",  o=1, dist="t"),        # asymmetry + fat tails
    "EGARCH-t":   fit(r, vol="EGARCH", o=1, dist="t"),        # log-variance asymmetry
    "GJR-skewt":  fit(r, vol="Garch",  o=1, dist="skewt"),    # + skew
}

for name, res in models.items():
    print(f"\n===== {name} =====")
    print(res.summary().tables[1])   # parameter table

Voor vol='EGARCH' regelt het argument o de asymmetrische ($\gamma z$) term en regelen p/q de magnitude- en lag-termen; o=1, p=1, q=1 is de standaard EGARCH(1,1). Eén valkuil: de parameternamen van EGARCH in arch zijn dezelfde letters, maar de tekenconventie op de asymmetrieterm is die van Nelson, dus een negatieve schatting is het hefboomeffect. We verifiëren dit vanuit de news impact curve in plaats van uit het geheugen.

De GJR-fit lezen

Een GJR-$t$-parametertabel ziet er ongeveer zo uit (illustratieve waarden, geen gerapporteerd experiment — herfit op je eigen data):

                  coef    std err       t      P>|t|
omega           0.0480     0.017     2.82     0.005
alpha[1]        0.0620     0.018     3.44     0.001
gamma[1]        0.0910     0.028     3.25     0.001
beta[1]         0.8850     0.021    42.1      0.000
nu              4.30       0.55      7.82     0.000

Hoe je het leest:

• gamma[1] = 0.091 met een $t$-statistiek boven 3 is een statistisch significant hefboomeffect. Na een negatieve schok is de impact van de gekwadrateerde schok $\alpha + \gamma = 0.062 + 0.091 = 0.153$; na een positieve schok is het slechts $\alpha = 0.062$. Slecht nieuws beweegt de volatiliteit van dit model ongeveer $2.5\times$ zo veel als goed nieuws van dezelfde grootte.
• nu = 4.3 bevestigt dikke staarten — ver van Gaussisch ($\nu \to \infty$), en laag genoeg dat het vierde moment nauwelijks eindig is. Een Gaussische VaR op deze reeks zou sterk te optimistisch zijn.
• De persistentie is $\alpha + \beta + \tfrac{1}{2}\gamma = 0.062 + 0.885 + 0.0455 \approx 0.993$ — zeer hoog, zoals gebruikelijk voor dagelijkse crypto: schokken dempen langzaam uit en volatiliteit is sterk geclusterd.

De belangrijkste regel om te controleren is de $\gamma$-rij. Als de $p$-waarde ervan groot is, verdient de asymmetrische term zijn plaats niet op dit activum en dit venster, en zou je het eenvoudigere symmetrische model moeten prefereren. Dit is modelselectiediscipline, geen versiersel — hieronder meer daarover.

Modellen vergelijken met informatiecriteria

De log-likelihood verbetert altijd wanneer je parameters toevoegt, dus je kunt niet alleen op log-likelihood selecteren. Gebruik AIC/BIC, die het aantal parameters bestraffen (BIC agressiever):

def compare(models: dict) -> pd.DataFrame:
    rows = []
    for name, res in models.items():
        rows.append({
            "model": name,
            "n_params": len(res.params),
            "loglik": res.loglikelihood,
            "AIC": res.aic,
            "BIC": res.bic,
        })
    df = pd.DataFrame(rows).set_index("model")
    return df.sort_values("BIC")

print(compare(models))

Vuistregels voor de interpretatie: een BIC-verbetering van meer dan ~6 ten opzichte van de baseline is sterk bewijs dat de extra structuur echt is; een verschil van 1–2 is ruis. Als GJR-t GARCH-N verslaat met 30+ BIC-punten, maar GJR-skewt GJR-t slechts met 1 verslaat, houd dan de $t$ aan en laat de scheefheid vallen — de scheefheidsparameter betaalt zichzelf niet terug op deze data. Lees AIC/BIC niet als vervanging voor out-of-sample-validatie; ze belonen in-sample-fit gecorrigeerd voor complexiteit, wat noodzakelijk maar niet voldoende is. De echte toets is de VaR-backtest en, uiteindelijk, walk-forward-evaluatie.

De news impact curve plotten

Dit is de beloningsplot — hij maakt asymmetrie zichtbaar en verifieert de EGARCH-tekenconventie.

import matplotlib.pyplot as plt

def news_impact_curve(res, shock_grid):
    """
    Next-period conditional variance as a function of the last shock,
    holding sigma_{t-1} at the model's unconditional level.
    Works for symmetric GARCH, GJR (o=1), and EGARCH.
    """
    p = res.params
    vol_name = res.model.volatility.__class__.__name__
    sigma2_bar = np.mean(res.conditional_volatility ** 2)

    omega = p["omega"]
    alpha = p.get("alpha[1]", 0.0)
    gamma = p.get("gamma[1]", 0.0)
    beta  = p.get("beta[1]", 0.0)

    if vol_name == "EGARCH":
        sig_prev = np.sqrt(sigma2_bar)
        z = shock_grid / sig_prev
        E_abs = np.sqrt(2 / np.pi)   # approx; arch uses the fitted dist's E|z|
        log_s2 = (omega + beta * np.log(sigma2_bar)
                  + alpha * (np.abs(z) - E_abs) + gamma * z)
        return np.exp(log_s2)
    else:
        ind = (shock_grid < 0).astype(float)
        return omega + (alpha + gamma * ind) * shock_grid**2 + beta * sigma2_bar

shocks = np.linspace(-8, 8, 401)   # daily % shocks
plt.figure(figsize=(8, 5))
for name in ["GARCH-N", "GJR-t", "EGARCH-t"]:
    nic = news_impact_curve(models[name], shocks)
    plt.plot(shocks, nic, label=name, lw=2)
plt.axvline(0, color="gray", lw=0.8, ls="--")
plt.xlabel("Last shock  $\\varepsilon_{t-1}$  (daily %)")
plt.ylabel("Next-period conditional variance  $\\sigma_t^2$")
plt.title("News impact curve: symmetric vs GJR vs EGARCH (BTC)")
plt.legend()
plt.tight_layout()

Wanneer je dit uitvoert, is de symmetrische GARCH-N-curve een zuivere parabool gecentreerd op nul — een schok van $-6\%$ en $+6\%$ geeft identieke variantie. GJR-t is een parabool met een knik in de oorsprong, hoger op de linkerarm. EGARCH-t is de exponentiële V, en als de linkerarm boven de rechterarm ligt, heb je in één oogopslag het hefboomeffect en de tekenconventie bevestigd. Als de linkerarm van EGARCH onder de rechter ligt, is ofwel $\gamma$ positief geschat (een up-vol-regime) of heb je het teken omgekeerd — de plot vertelt je welke van de twee zonder enig raden.

Een zij-aan-zij-vergelijking van de vier modellen

Voordat we ons tot risico wenden, helpt het om de vier modellen naast elkaar te houden. Elke rij is een ontwerpbeslissing, en de kolommen laten zien wat die beslissing kost en oplevert.

Het patroon om te internaliseren: de sprong van kolom 1 naar kolom 2 — het tegelijk toevoegen van zowel asymmetrie als dikke staarten — is waar bijna alle verbetering in risicokalibratie zich bevindt. De daaropvolgende verfijningen (de functionele vorm van EGARCH, de scheefheidsterm) zijn echt maar tweederangs, en op veel crypto-reeksen vallen ze binnen de ruis. Besteed je modelbudget aan de eerste sprong en wees sceptisch over de rest.

Risicotoepassing: VaR en Expected Shortfall

Een chiquer volatiliteitsmodel fitten is alleen de moeite waard als het een beslissing verbetert. De schoonste beslissing om te verbeteren is de eenstaps-staartrisicoforecast: hoe slecht kan morgen zijn? We produceren een eendaags-vooruit Value-at-Risk en Expected Shortfall (ook wel Conditional VaR, die de HRP/CVaR-portefeuillepijplijn als doelfunctie gebruikt) rechtstreeks uit de gefitte GARCH-$t$/skew-$t$-forecast.

Van conditionele verdeling naar VaR

De GARCH-machinerie geeft een eenstaps-forecast van het conditionele gemiddelde $\mu_{t+1}$ en de conditionele standaarddeviatie $\sigma_{t+1}$. Het rendement wordt gemodelleerd als $r_{t+1} = \mu_{t+1} + \sigma_{t+1} z_{t+1}$ met $z_{t+1}$ getrokken uit de gefitte gestandaardiseerde verdeling (Gaussisch, $t$ of skew-$t$). Dus het $\alpha$-kwantiel van het rendement is slechts een affiene transformatie van het $\alpha$-kwantiel van de gestandaardiseerde verdeling:

$\text{VaR}_{\alpha}(t+1) = -\Bigl(\mu_{t+1} + \sigma_{t+1}\, F_z^{-1}(1-\alpha)\Bigr)$

waarbij $F_z^{-1}$ het kwantiel (inverse CDF) van de gestandaardiseerde innovatie is en het leidende minteken de conventie volgt dat VaR een positief verliesgetal is. Voor een 99%-VaR is $\alpha = 0.99$ en vul je $F_z^{-1}(0.01)$ in. Het hele voordeel van de $t$/skew-$t$ komt hier tot uiting: $F_z^{-1}(0.01)$ is negatiever dan de Gaussische $-2.326$, dus de VaR is eerlijk groter.

Expected Shortfall

VaR vertelt je de drempel; het zegt niets over hoe slecht de doorbraak is wanneer die zich voordoet. Expected Shortfall — het gemiddelde verlies conditioneel op het overschrijden van VaR — doet dat wel, en het is coherent (subadditief), wat de reden is dat het de risicomaat is achter CVaR-optimalisatie en waarom Basel ernaar is overgestapt. Voor een locatie-schaalmodel,

$\text{ES}_{\alpha}(t+1) = -\Bigl(\mu_{t+1} + \sigma_{t+1}\, \mathbb{E}\!\left[z \mid z \leq F_z^{-1}(1-\alpha)\right]\Bigr).$

De conditionele-staartverwachtingsterm $\mathbb{E}[z \mid z \le q]$ heeft gesloten vormen voor de standaardverdelingen. Voor de Gaussische, met $q = \Phi^{-1}(1-\alpha)$,

$\mathbb{E}[z \mid z \le q] = -\frac{\phi(q)}{1-\alpha},$

waarbij $\phi$ de standaardnormale pdf is. Voor de gestandaardiseerde Student-$t$ met $\nu$ vrijheidsgraden en $q = t_\nu^{-1}(1-\alpha)$ (op de gestandaardiseerde schaal), is de staartverwachting

$\mathbb{E}[z \mid z \le q] = -\frac{\nu + q^2}{\nu - 1}\cdot\frac{g_\nu(q)}{1-\alpha},$

waarbij $g_\nu$ de gestandaardiseerde-$t$ pdf is. De $t$-Expected-Shortfall overtreft de Gaussische met meer dan de VaR dat doet, omdat de $t$-staart niet alleen verder naar buiten ligt — hij is dikker, zodat het gemiddelde verlies voorbij de drempel onevenredig groot is. Die extra kloof is het getal dat een Gaussisch model voor je verbergt.

VaR en ES berekenen uit een gefit arch-model

arch-verdelingen stellen een ppf (kwantiel)-methode beschikbaar, zodat we het gestandaardiseerde kwantiel direct kunnen verkrijgen en niets opnieuw hoeven te herleiden. Voor ES integreren we numeriek, wat robuust is en uniform werkt over normal/t/skewt.

from scipy import integrate

def var_es_forecast(res, alpha=0.99):
    """
    One-step-ahead VaR and ES at level alpha, on the same (x100) scale
    as the returns fed to the model. Divide by 100 for fractional units.
    """
    fc = res.forecast(horizon=1, reindex=False)
    mu = fc.mean.iloc[-1, 0]
    sigma = np.sqrt(fc.variance.iloc[-1, 0])

    dist = res.model.distribution          # StudentsT / SkewStudent / Normal
    dp = [res.params[k] for k in res.params.index
          if k in ("nu", "eta", "lambda", "lam")]

    z_q = dist.ppf(1 - alpha, dp) if dp else dist.ppf(1 - alpha)
    z_q = float(np.atleast_1d(z_q)[0])

    def pdf(z):
        arr = np.atleast_1d(z).astype(float)
        lp = dist.loglikelihood(dp, arr, np.ones_like(arr), individual=True) \
             if dp else dist.loglikelihood([], arr, np.ones_like(arr),
                                           individual=True)
        return np.exp(lp)

    num, _ = integrate.quad(lambda z: z * pdf(z), -30, z_q, limit=200)
    es_z = num / (1 - alpha)               # E[z | z <= z_q]

    var = -(mu + sigma * z_q)
    es  = -(mu + sigma * es_z)
    return {"mu": mu, "sigma": sigma, "z_q": z_q,
            "VaR": var, "ES": es}

for name in ["GARCH-N", "GJR-t", "GJR-skewt"]:
    out = var_es_forecast(models[name], alpha=0.99)
    print(f"{name:11s}  sigma={out['sigma']:.2f}%  "
          f"z_q={out['z_q']:+.2f}  "
          f"VaR99={out['VaR']:.2f}%  ES99={out['ES']:.2f}%")

De z_q-kolom is het hele verhaal in één getal. Het Gaussische model gebruikt $z_q \approx -2.33$; de $t$ met $\nu \approx 4.3$ gebruikt iets in de buurt van $-3.3$; de skew-$t$ duwt het linker-kwantiel nog verder naar buiten terwijl het rechter naar binnen wordt getrokken. Dezelfde $\sigma_{t+1}$, wezenlijk grotere VaR. Als je Gaussische VaR op crypto hebt gedraaid, is dit de kloof die je stilletjes hebt geabsorbeerd.

Eén stap versus meerdere stappen: een kanttekening

Alles hierboven is een eendaags-vooruit forecast, en dat is waar GARCH-VaR het schoonst is. Twee zaken compliceren langere horizonnen en je zou ze moeten kennen voordat je extrapoleert.

Ten eerste keren variantieforecasts terug naar het gemiddelde. De $h$-staps-vooruit conditionele variantie uit een stationaire GARCH convergeert naar het onconditionele niveau $\bar\sigma^2$ naarmate $h$ groeit, en de cumulatieve $h$-dagse variantie is de som van de per-stap-forecasts — het is niet $h\times\sigma_{t+1}^2$ tenzij de volatiliteit op zijn langetermijngemiddelde is. De naïeve "wortel-van-de-tijd"-schaling $\text{VaR}(h) = \sqrt{h}\,\text{VaR}(1)$ negeert deze terugkeer naar het gemiddelde en is precies fout na een schok, wanneer je het getal het hardst nodig hebt. Gebruik het eigen meerstaps-variantiepad van het model.

Ten tweede is de verdeling van een meerdaags rendement niet dezelfde vorm als de eendaagse innovatie. Het optellen van meerdere $t$-verdeelde dagelijkse schokken (via de niet-lineaire GARCH-recursie) geeft geen $t$-verdeling op de $h$-daagse horizon; er is geen schone gesloten vorm. Voor meerdaagse VaR is de eerlijke route simulatie: trek innovatiepaden uit de gefitte gestandaardiseerde verdeling, laat ze door de GARCH-recursie lopen om gesimuleerde rendementspaden te verkrijgen, aggregeer tot $h$-daagse rendementen, en lees het empirische kwantiel af. Dat handelt ook op natuurlijke wijze het skew-$t$-geval af, waar helemaal geen analytisch meerhorizon-kwantiel bestaat. De eenstaps-analytische formules in deze post zijn exact; behandel elke meerstaps-kortere weg als een benadering die je moet valideren.

VaR backtesten: Kupiec en Christoffersen

Een VaR-forecast is een probabilistische claim: "het verlies zal deze drempel op slechts $(1-\alpha)$ van de dagen overschrijden." Je toetst het door overtredingen te tellen (dagen waarop het gerealiseerde verlies de forecast-VaR overschreed) over een walk-forward-evaluatie en twee dingen te controleren. Ten eerste: is de overtredingsgraad juist? Ten tweede: zijn de overtredingen onafhankelijk, of clusteren ze (wat betekent dat het model precies faalt wanneer het ertoe doet, tijdens volatiliteitspieken)?

Laat $I_t = \mathbf{1}\{\text{loss}_t > \text{VaR}_t\}$ de overtredingsreeks zijn, $N = \sum I_t$ het aantal overtredingen over $T$ dagen, en $\hat{\pi} = N/T$ de waargenomen graad. Doelgraad $p = 1-\alpha$.

Kupiecs onconditionele-dekkingstoets (1995) controleert $\hat\pi \approx p$ via een likelihood-ratio:

$LR_{uc} = -2\log\!\left[\frac{p^{N}(1-p)^{T-N}}{\hat\pi^{N}(1-\hat\pi)^{T-N}}\right] \;\sim\; \chi^2_1.$

Christoffersens onafhankelijkheidstoets (1998) controleert dat een overtreding vandaag niet wordt voorspeld door een overtreding gisteren. Laat $n_{ij}$ de overgangen van toestand $i$ naar toestand $j$ in de overtredingsreeks tellen, $\pi_{01} = n_{01}/(n_{00}+n_{01})$, $\pi_{11} = n_{11}/(n_{10}+n_{11})$, en $\pi = (n_{01}+n_{11})/T$. Dan

$LR_{ind} = -2\log\!\left[\frac{(1-\pi)^{n_{00}+n_{10}}\,\pi^{\,n_{01}+n_{11}}}{(1-\pi_{01})^{n_{00}}\pi_{01}^{n_{01}}(1-\pi_{11})^{n_{10}}\pi_{11}^{n_{11}}}\right] \;\sim\; \chi^2_1.$

De twee combineren tot de conditionele-dekkingstoets $LR_{cc} = LR_{uc} + LR_{ind} \sim \chi^2_2$, die tegelijk correcte graad én onafhankelijkheid controleert. Een model kan Kupiec doorstaan (juiste aantal overtredingen) en toch Christoffersen falen (ze hoopten zich allemaal op in één crashweek) — dat is de faalmodus die je het liefst wilt betrappen, omdat geclusterde overtredingen degene zijn die een rekening opblazen.

from scipy.stats import chi2

def var_backtest(losses, var, p):
    """
    losses, var: 1D arrays, same length; loss>0 is a loss, var>0 threshold.
    p = 1 - alpha (target violation rate, e.g. 0.01 for 99% VaR).
    Returns Kupiec, Christoffersen-independence, and conditional-coverage LR.
    """
    losses = np.asarray(losses)
    var = np.asarray(var)
    I = (losses > var).astype(int)          # violation indicators
    T = len(I)
    N = int(I.sum())
    pi_hat = N / T

    eps = 1e-12
    ll_null = N * np.log(p + eps) + (T - N) * np.log(1 - p + eps)
    ll_alt  = N * np.log(pi_hat + eps) + (T - N) * np.log(1 - pi_hat + eps)
    LR_uc = -2 * (ll_null - ll_alt)
    p_uc = 1 - chi2.cdf(LR_uc, df=1)

    n00 = np.sum((I[:-1] == 0) & (I[1:] == 0))
    n01 = np.sum((I[:-1] == 0) & (I[1:] == 1))
    n10 = np.sum((I[:-1] == 1) & (I[1:] == 0))
    n11 = np.sum((I[:-1] == 1) & (I[1:] == 1))
    pi01 = n01 / (n00 + n01) if (n00 + n01) else 0.0
    pi11 = n11 / (n10 + n11) if (n10 + n11) else 0.0
    pi   = (n01 + n11) / (n00 + n01 + n10 + n11)

    def _ll(pi01, pi11, pi_pooled, use_pooled):
        if use_pooled:
            a = (n00 + n10) * np.log(1 - pi_pooled + eps)
            b = (n01 + n11) * np.log(pi_pooled + eps)
            return a + b
        return (n00 * np.log(1 - pi01 + eps) + n01 * np.log(pi01 + eps)
                + n10 * np.log(1 - pi11 + eps) + n11 * np.log(pi11 + eps))

    LR_ind = -2 * (_ll(pi01, pi11, pi, True) - _ll(pi01, pi11, pi, False))
    p_ind = 1 - chi2.cdf(LR_ind, df=1)

    LR_cc = LR_uc + LR_ind
    p_cc = 1 - chi2.cdf(LR_cc, df=2)

    return {
        "T": T, "violations": N,
        "obs_rate": pi_hat, "target_rate": p,
        "LR_uc": LR_uc, "p_uc": p_uc,
        "LR_ind": LR_ind, "p_ind": p_ind,
        "LR_cc": LR_cc, "p_cc": p_cc,
    }

Om de losses/var-inputs eerlijk te genereren, herfit je (of forecast op zijn minst opnieuw) op een uitdijend of rollend venster en registreer je de eenstaps-vooruit VaR voor elke out-of-sample-dag, en vergelijk je die vervolgens met het gerealiseerde verlies voor die dag. Backtest VaR nooit in-sample — een model dat gefit is op dezelfde crash die het geacht wordt te voorspellen, ziet er veel beter uit dan het is. Dit is dezelfde discipline als backtest-live-pariteit: de evaluatie mag alleen informatie gebruiken die beschikbaar is op het beslismoment.

def walk_forward_var(returns, alpha=0.99, start=750, refit_every=25,
                     vol="Garch", o=1, dist="t"):
    """
    Expanding-window one-step VaR. Refit every `refit_every` days
    (refitting daily is correct but slow; every ~25 days is a common
    compromise -- validate the shortcut on your data).
    """
    losses, vars_ = [], []
    res = None
    for t in range(start, len(returns)):
        if res is None or (t - start) % refit_every == 0:
            am = arch_model(returns.iloc[:t], mean="Constant",
                            vol=vol, p=1, o=o, q=1, dist=dist)
            res = am.fit(disp="off")
        fc = res.forecast(horizon=1, reindex=False)
        mu = fc.mean.iloc[-1, 0]
        sig = np.sqrt(fc.variance.iloc[-1, 0])
        dp = [res.params[k] for k in res.params.index
              if k in ("nu", "eta", "lambda", "lam")]
        z_q = float(np.atleast_1d(
            res.model.distribution.ppf(1 - alpha, dp) if dp
            else res.model.distribution.ppf(1 - alpha))[0])
        var_t = -(mu + sig * z_q)
        vars_.append(var_t)
        losses.append(-returns.iloc[t])     # realized loss for that day
    return np.array(losses), np.array(vars_)

losses, vars_ = walk_forward_var(r, alpha=0.99, dist="t")
print(var_backtest(losses, vars_, p=0.01))

De lezing: een goed gekalibreerde 99%-VaR toont een waargenomen graad in de buurt van 1%, een niet-significante Kupiec (grote p_uc), en een niet-significante Christoffersen (grote p_ind) — geen clustering. In de praktijk is het eerlijke resultaat op crypto dat GARCH-Normal Kupiec faalt (te veel overtredingen, p_uc piepklein) terwijl GJR-$t$ of EGARCH-$t$ slaagt of dicht in de buurt komt. Dat contrast is het hele argument van deze post weergegeven als een hypothesetoets. Als zelfs het $t$-model geclusterde overtredingen vertoont (kleine p_ind), is je volatiliteitsdynamiek nog steeds verkeerd gespecificeerd — vaak een teken dat je een langer geheugen (component/FIGARCH) of een regimelaag nodig hebt, wat aansluit bij regimedetectie met HMM's.

Modellen rangschikken op staartverlies, niet alleen slagen/falen

Kupiec en Christoffersen geven je een binair oordeel — het model wordt wel of niet verworpen. Dat is noodzakelijk maar grof: twee modellen kunnen beide "slagen" terwijl het ene wezenlijk scherper is. Om concurrerende VaR-forecasts te rangschikken, scoor je ze met een strikt consistente verliesfunctie voor het kwantiel, de pinball- (kwantiel-)loss:

$L_\tau(r, q) = \bigl(\tau - \mathbf{1}\{r < q\}\bigr)\,(r - q), \qquad \tau = 1-\alpha,$

waarbij $q$ het (met teken) VaR-kwantiel is en $r$ het gerealiseerde rendement. Gemiddeld over de out-of-sample-dagen betekent een lagere gemiddelde pinball-loss een beter gekalibreerd en scherper kwantiel; omdat de loss consistent is voor het $\tau$-kwantiel, beloont het minimaliseren ervan een model niet voor het lui-breed zijn. Om twee modellen formeel te vergelijken, voer je hun per-dag-verliesverschillen in een Diebold-Mariano-toets.

def pinball_loss(returns, var, alpha=0.99):
    tau = 1 - alpha
    q = -np.asarray(var)               # VaR is a positive loss; quantile is negative
    r = np.asarray(returns)
    hit = (r < q).astype(float)
    return np.mean((tau - hit) * (r - q))

Voor Expected Shortfall in het bijzonder: merk op dat ES niet op zichzelf eliciteerbaar is (er bestaat geen verliesfunctie waarvan de minimalisator ES alleen is), wat een echte theoretische rimpel is: je evalueert ES gezamenlijk met VaR met behulp van de Fissler-Ziggel-scoringsregels, of je valt terug op de eenvoudigere praktijk van controleren dat de gemiddelde doorbraakgrootte overeenkomt met de door het model voorspelde ES. Een grove maar bruikbare ES-controle: vergelijk onder de VaR-overtredingsdagen het gemiddelde gerealiseerde verlies met de gemiddelde forecast-ES op die dagen — ze zouden dicht bij elkaar moeten liggen.

De regelgevende kadering is de Basel-verkeerslicht-aanpak: over 250 handelsdagen is 0-4 overtredingen van een 99%-VaR "groen" (aanvaardbaar), 5-9 is "geel" (extra aandacht), 10+ is "rood" (het model wordt verworpen en de kapitaalmultiplicatoren stijgen). Het is een grovere neef van Kupiec, maar het is de taal die risicocommissies daadwerkelijk spreken, en het is de moeite waard om naast de LR-statistieken te rapporteren.

Praktische overwegingen

Wanneer de extra parameters zich niet uitbetalen

De eerlijke standaard is scepsis ten opzichte van complexiteit. Elke parameter die je toevoegt is een knop die de optimizer kan overfitten, en asymmetrische dikstaartige GARCH heeft er meerdere. Concrete richtlijnen:

• Illiquide of korte steekproeven. Met een paar honderd dagelijkse observaties zal de standaardfout op $\gamma$ en $\lambda$ groot zijn, en zul je asymmetrieën "detecteren" die steekproefruis zijn. Op een nieuwe of dunne altcoin is een symmetrische GARCH-$t$ vaak het meest complexe model dat de data kan ondersteunen. Skew-$t$-EGARCH fitten op 200 dagen is jezelf voor de gek houden.
• De scheefheidsterm verdient zijn kosten vaak niet. In de praktijk is de overgang van Normal → $t$ een grote, betrouwbare verbetering (dikke staarten zijn echt en sterk). De overgang van $t$ → skew-$t$ is vaak marginaal — een BIC-winst van 1 of 2, soms negatief. Voeg scheefheid alleen toe wanneer de data er duidelijk om vragen.
• EGARCH versus GJR is meestal een gelijkspel op dagelijkse data. Ze coderen hetzelfde kwalitatieve verhaal met verschillende functionele vormen. Kies op basis van de out-of-sample VaR-backtest, niet op basis van welke de mooiere in-sample log-likelihood heeft.
• Hogere frequentie verschuift het antwoord. Op uur- of minuutbalken domineren intraday-seizoensinvloeden en microstructuur, en een gewone GARCH in dagstijl is verkeerd gespecificeerd ongeacht asymmetrie. Ander probleem, ander gereedschap.

Dit is dezelfde les als eerlijke evaluatie zonder robuuste edge: een complexer model dat out-of-sample-toetsing niet overleeft, is slechter dan het eenvoudige model dat het verving, omdat het de illusie van precisie met zich meedraagt. Rapporteer het negatieve resultaat — "scheefheid hielp niet op ETH" — als een echte bevinding, en gebruik walk-forward-optimalisatie als scheidsrechter, niet in-sample-AIC.

Dit zijn de marginalen waarop iedereen verder bouwt

De modellen hier zijn geen eindpunt; ze zijn de univariate bouwsteen voor de machinerie van gezamenlijk risico. De post copula-modellen voor gezamenlijk crypto-risico gebruikt precies EGARCH/GJR-$t$ als de GARCH-EVT-marginalen voordat een vine-copula wordt gefit — je fit een asymmetrische dikstaartige GARCH per activum, extraheert gestandaardiseerde residuen, en modelleert pas dan de cross-asset-afhankelijkheid. Als je marginaal een symmetrische Gaussische GARCH is, erft de copula de staartfouten ervan, hoe goed het afhankelijkheidsmodel ook is. Rommelmarginalen, rommel-gezamenlijke-VaR.

Voor het multivariate volatiliteitsprobleem — tijdvariërende correlaties in plaats van per-activum-variantie — zie Deel 3, DCC-GARCH, dat een dynamisch-correlatiemodel bovenop deze univariate fits legt. En voor het omzetten van een volatiliteitsforecast in positiegrootte en een trading-backtest, gebruikt Deel 4 over volatiliteitstargeting de $\sigma_{t+1}$-forecasts uit precies deze modellen om de blootstelling omgekeerd te schalen ten opzichte van het voorspelde risico.

Een verdelingsvrij alternatief

Alles in de risicosectie rust op een parametrische aanname: dat gestandaardiseerde residuen een $t$ of skew-$t$ volgen. Die aanname is toetsbaar en meestal redelijk, maar hij kan falen. Als je liever helemaal geen staartvorm wilt vastleggen, geeft conformal prediction verdelingsvrije voorspellingsintervallen met eindige-steekproef-dekkingsgaranties — een werkelijk andere filosofie die geen enkele claim doet over de innovatieverdeling. De twee benaderingen zijn complementair: parametrische GARCH-$t$ geeft je een volledige conditionele dichtheid (en dus ES, wat conformale intervallen niet rechtstreeks bieden), terwijl conformal je dekking geeft die standhoudt zelfs wanneer je dichtheid verkeerd is. In productie is beide gebruiken als kruiscontrole goedkope verzekering.

Numerieke en workflow-hygiëne

• Schaal rendementen met 100. GARCH-optimizers convergeren veel betrouwbaarder op procentrendementen dan op ruwe fractionele rendementen. Vergeet niet VaR/ES te ontschalen als je in fractionele eenheden rapporteert.
• Let op de persistentie. Als $\alpha + \beta + \tfrac12\gamma$ boven ~0.999 schat, is het model bijna-geïntegreerd (IGARCH-achtig); forecasts keren extreem langzaam terug naar het gemiddelde en langehorizon-variantieforecasts worden onbetrouwbaar. Niet noodzakelijk verkeerd, maar signaleer het.
• Convergentiefouten op rollende vensters. De logvorm van EGARCH vermijdt positiviteitsbeperkingen, maar kan nog steeds falen te convergeren op een pathologisch venster. Verpak fit() in een try/except en val terug op de parameters van het vorige venster in plaats van een live backtest te laten crashen.
• Gemiddelde-model. We hebben overal een constant gemiddelde gebruikt. Voor de meeste dagelijkse crypto is het conditionele gemiddelde dicht bij nul en wordt het overspoeld door volatiliteit; besteed geen modelcomplexiteit aan het proberen te voorspellen ervan tenzij je een echte reden hebt.

Samenvatting

• Gewone GARCH(1,1) heeft twee structurele gebreken: het is symmetrisch (reageert identiek op $+x\%$ en $-x\%$ omdat schokken als $\varepsilon^2$ binnenkomen) en het gaat uit van Gaussische innovaties (waardoor het de dikke staarten van crypto onderschat). Beide kosten echt geld via optimistische VaR.
• GJR-GARCH voegt een drempelterm $\gamma\,I(\varepsilon_{t-1}<0)\,\varepsilon_{t-1}^2$ toe. Een significante $\gamma > 0$ is het hefboomeffect: slecht nieuws verhoogt de volatiliteit meer. Positiviteit vereist $\alpha+\gamma\ge0$; persistentie is $\alpha+\beta+\tfrac12\gamma$.
• EGARCH modelleert $\log\sigma_t^2$, dus geen positiviteitsbeperkingen en stationariteit is slechts $|\beta|<1$. Asymmetrie komt binnen via een term met teken $\gamma z_{t-1}$ (hefboom is $\gamma<0$ in deze conventie) gescheiden van een magnitude-term $|z_{t-1}|$.
• De news impact curve — variantie van de volgende periode tegen de laatste schok — maakt de asymmetrie zichtbaar en verifieert de EGARCH-tekenconventie in één oogopslag.
• Student-$t$-innovaties (dist='t') corrigeren de staarten via een vrijheidsgraden $\nu$ (typisch 3–6 voor crypto); Hansens skew-$t$ (dist='skewt') voegt een scheefheid $\lambda$ toe voor een zwaardere linkerstaart. De overgang Normal → $t$ is een grote betrouwbare winst; $t$ → skew-$t$ is vaak marginaal.
• VaR en ES volgen uit de gefitte conditionele verdeling: $\text{VaR}_\alpha = -(\mu_{t+1} + \sigma_{t+1}F_z^{-1}(1-\alpha))$, waarbij het dikstaartige kwantiel het risico eerlijk groter maakt dan Gaussisch. ES (coherent, $\approx$ CVaR) vangt het gemiddelde verlies voorbij VaR.
• Backtest met Kupiec en Christoffersen. Kupiec controleert de overtredingsgraad; Christoffersen controleert dat overtredingen niet geclusterd zijn. Een model kan het ene doorstaan en het andere falen — geclusterde overtredingen zijn de gevaarlijke faalmodus. Backtest strikt out-of-sample.
• Discipline boven complexiteit. Voeg asymmetrie/scheefheid alleen toe wanneer het BIC én een out-of-sample VaR-backtest overleeft. Op korte of illiquide reeksen wint het eenvoudigere model meestal.

References:

• Nelson, D. B. (1991). Conditional heteroskedasticity in asset returns: A new approach. Econometrica, 59(2), 347-370. DOI
• Glosten, L. R., Jagannathan, R., & 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. DOI
• Bollerslev, T. (1987). A conditionally heteroskedastic time series model for speculative prices and rates of return. Review of Economics and Statistics, 69(3), 542-547. DOI
• Hansen, B. E. (1994). Autoregressive conditional density estimation. International Economic Review, 35(3), 705-730. DOI
• Engle, R. F., & Ng, V. K. (1993). Measuring and testing the impact of news on volatility. Journal of Finance, 48(5), 1749-1778. DOI
• Christoffersen, P. F. (1998). Evaluating interval forecasts. International Economic Review, 39(4), 841-862. DOI
• Kupiec, P. H. (1995). Techniques for verifying the accuracy of risk measurement models. Journal of Derivatives, 3(2), 73-84. DOI
• Black, F. (1976). Studies of stock price volatility changes. Proceedings of the American Statistical Association, Business and Economic Statistics Section, 177-181.
• McNeil, A. J., Frey, R., & Embrechts, P. (2015). Quantitative Risk Management: Concepts, Techniques and Tools (2nd ed.). Princeton University Press.
• Sheppard, K. (2023). arch: Autoregressive conditional heteroskedasticity models in Python. GitHub.