GARCH models in Python: ARCH, GARCH, GJR and EGARCH

Asset returns are hard to predict, but their volatility is not. Calm days tend to follow calm days and turbulent days follow turbulent ones, a pattern called volatility clustering. A regression with a constant error variance ignores this: its forecast intervals are too wide in quiet periods and too narrow in crises, and it says nothing about risk, which is often the quantity of interest. ARCH-family models let the variance of today's shock depend on the size of recent shocks and on recent variances.

What GARCH models estimate

The mean equation is an ordinary regression, often just a constant. Its error is the product of a conditional standard deviation and an independent shock with mean 0 and variance 1. The variance equation describes how the conditional variance h_t evolves. In the GARCH(1,1) model:

y_t = x_t'b + e_t,        e_t = √h_t · z_t,        z_t ~ iid (0, 1)
h_t = ω + α e²_{t−1} + β h_{t−1}

Engle (1982) introduced the ARCH model, in which h_t depends on past squared shocks only. Bollerslev (1986) added lagged variances, the GARCH term, which captures long-lasting volatility with few parameters. The sum α + β is the persistence: the closer it is to one, the more slowly a volatility shock fades. When it is below one, the variance reverts to the unconditional level ω / (1 − α − β).

The main variants change the variance equation or the shock distribution:

• GJR-GARCH (Glosten, Jagannathan and Runkle 1993) adds a term γ e²_{t−1} for negative shocks only, so bad news can raise volatility more than good news of the same size. This leverage effect is typical of equity returns.
• EGARCH (Nelson 1991) models the logarithm of the variance with separate terms for the sign and the size of the standardized shock. The variance is positive for any parameter values.
• Power ARCH, IGARCH and ARCH-in-mean respectively estimate the power of the variance equation, impose a persistence of exactly one, or let the conditional variance enter the mean equation as a risk premium.
• Student-t or GED shocks allow fatter tails than the normal distribution.

All are estimated by maximum likelihood. A sensible workflow is to test the residuals of the mean equation for ARCH effects, start with GARCH(1,1), add an asymmetric term or fat-tailed shocks when the data call for them, and compare candidates by information criteria and by diagnostics of the standardized residuals z_t = e_t / √h_t.

In OpenEconometrics

oe.arch fits the whole family. The call shape is oe.arch(data=df, y="ret", time="day", model="gjr"). With the data declared by tsset day, the GJR model in the example below corresponds to the Stata command arch ret, arch(1) tarch(1) garch(1), and the default GARCH(1,1) to arch ret, arch(1) garch(1).

• Time. time= is optional. Rows are sorted by it and integer periods must be consecutive. A datetime column is ranked and its observations are treated as consecutive, so weekends and holidays of daily data are simply skipped; a warning records this. Without time=, the input order is the time order.
• Variance equation. model= is "arch", "garch" (the default), "gjr" (alias "tarch"), "egarch", "parch", "igarch", "cgarch" or "figarch". arch= and garch= take an order (arch=2 means lags 1 and 2) or a list of lags (arch=[1, 3]). At least one ARCH term is required.
• Shocks. dist= is "normal", "t", "ged", "skewt" or "sged". The Student-t degrees of freedom are estimated as /lndfm2 = ln(ν − 2) and reported on their natural scale in fit.extra["distribution"].
• Mean equation. x= adds regressors, ar= and ma= add ARMA disturbances and archm="variance", "sd" or "log" adds an ARCH-in-mean term. variance_x= puts regressors in the variance equation, like Stata's het().
• Covariance. covariance="opg" (the default), "nonrobust" (observed information) or "robust", the quasi-maximum-likelihood sandwich, which remains valid when the shock distribution is misspecified. Tests are z tests.

Two conventions matter when comparing coefficients. The threshold term ARCH:L1.tarch multiplies the indicator of a negative shock, as in the original GJR paper. Stata's tarch() term multiplies the indicator of a positive shock, so Stata's arch coefficient corresponds to arch + tarch here and its tarch coefficient to −tarch; the fitted variances are the same. Second, no positivity or stationarity constraints are imposed on the variance coefficients. A persistence of one or more and negative variance coefficients are reported with a warning.

The result reports persistence and unconditional_variance in fit.metrics, and four diagnostics of the standardized residuals in fit.tests: Ljung–Box tests of z_t and z_t², an ARCH-LM test and a Jarque–Bera test. oe.forecast(fit, steps=h) (alias oe.arch_forecast) continues the fitted recursions and returns the mean forecast, the variance forecast, its square root and an interval. For GARCH and GJR the variance forecast is the exact conditional expectation; it approaches the unconditional variance at the rate of the persistence. To test a series for ARCH effects before modelling, use oe.archlm(df, "resid", lags=4). The ARCH family guide gives the likelihoods, presample values and covariance formulas.

Supported options

Options of oe.arch

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 2,500 daily returns with mean 0.04 and normal GJR-GARCH(1,1) shocks: ω = 0.02, α = 0.03, a threshold coefficient γ = 0.10 on negative shocks and β = 0.88. The persistence is 0.03 + 0.10/2 + 0.88 = 0.96, because half the shocks are negative, and the unconditional variance is 0.02 / 0.04 = 0.5. The script fits a symmetric GARCH(1,1) and the GJR model, compares them by AIC and forecasts five days ahead.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=2500, seed=31):
    """Daily returns with GJR-GARCH(1,1) variance: bad news raises volatility more than good news."""
    rng = np.random.default_rng(seed)
    z = rng.normal(size=n)
    e, h = np.zeros(n), np.zeros(n)
    h_prev, e_prev = 0.5, 0.0                     # start at the unconditional variance 0.02 / 0.04
    for t in range(n):
        h[t] = 0.02 + 0.03 * e_prev**2 + 0.10 * e_prev**2 * (e_prev < 0) + 0.88 * h_prev
        e[t] = np.sqrt(h[t]) * z[t]
        h_prev, e_prev = h[t], e[t]
    return pd.DataFrame({"day": np.arange(1, n + 1), "ret": 0.04 + e})


if __name__ == "__main__":
    df = simulate()
    garch = oe.arch(data=df, y="ret", time="day")               # GARCH(1,1), the default
    gjr = oe.arch(data=df, y="ret", time="day", model="gjr")    # adds the threshold term
    print(f"AIC  GARCH(1,1): {garch.metrics['aic']:.2f}  GJR-GARCH(1,1): {gjr.metrics['aic']:.2f}")
    print(gjr.summary())
    print(oe.forecast(gjr, steps=5).round(4))

With openecon 0.3.20a1 from PyPI the script prints:

AIC  GARCH(1,1): 5424.68  GJR-GARCH(1,1): 5401.59
GJR-GARCH(1,1) regression, normal errors — ret
Observations: 2500  |  Covariance: opg  |  Confidence: 95%

Term             Estimate  Std. error        z     P>|stat|    CI lower   CI upper
--------------  ---------  ----------  -------  -----------  ----------  ---------
[ret]
Intercept       0.0352286   0.0133843  2.63208   0.00848636  0.00899583  0.0614614
[ARCH]
ARCH:L1.arch    0.0435226   0.0146175  2.97744   0.00290663    0.014873  0.0721723
ARCH:L1.tarch    0.102992   0.0205339  5.01573  5.28331e-07   0.0627466   0.143238
ARCH:L1.garch    0.869741   0.0189233  45.9614            0    0.832652    0.90683
ARCH:Intercept  0.0204226  0.00506261  4.03401  5.48331e-05   0.0105001  0.0303451

log_likelihood: -2695.79  |  aic: 5401.59  |  bic: 5430.71  |  persistence: 0.96476
unconditional_variance: 0.57953  |  iterations: 11
Ljung-Box Q(40) test of the standardized residuals: chi2(40) = 38.365, p = 0.544
Ljung-Box Q(40) test of the squared standardized residuals: chi2(40) = 48.3895, p = 0.1703
ARCH-LM(5) test of the standardized residuals: chi2(5) = 4.91084, p = 0.4269
Jarque-Bera normality test of the standardized residuals: chi2(2) = 1.37218, p = 0.5035
   period  mean_forecast  variance_forecast  std_error  ci_low  ci_high
0    2501         0.0352             0.4086     0.6392 -1.2176   1.2881
1    2502         0.0352             0.4146     0.6439 -1.2268   1.2973
2    2503         0.0352             0.4204     0.6484 -1.2356   1.3061
3    2504         0.0352             0.4261     0.6527 -1.2441   1.3145
4    2505         0.0352             0.4315     0.6569 -1.2522   1.3226

Reading the output

• The asymmetric model wins. AIC falls from 5424.68 for GARCH(1,1) to 5401.59 for GJR, a gain of 23 points for one extra parameter.
• The leverage effect is recovered. The threshold coefficient is 0.103 (standard error 0.021, z = 5.02), against a true 0.10. A negative shock raises tomorrow's variance by 0.044 + 0.103 = 0.147 per unit of squared shock, a positive one by 0.044 only. The other estimates, 0.0435 for arch, 0.870 for garch and 0.0204 for the constant, all have intervals that contain the true 0.03, 0.88 and 0.02.
• Volatility is persistent. The estimated persistence is 0.965 (true 0.96). A volatility shock therefore loses half its effect in about ln 0.5 / ln 0.965 ≈ 19 trading days.
• The unconditional variance is imprecise. It is 0.580 against a true 0.5. The ratio ω / (1 − persistence) is sensitive to small errors in a persistence close to one.
• The model captures the clustering. The Ljung–Box test of the squared standardized residuals (p = 0.17) and the ARCH-LM test (p = 0.43) find no remaining ARCH effects, and Jarque–Bera (p = 0.50) is consistent with the normal shocks of the simulation.
• The forecasts revert to the mean. The sample ends in a calm spell, so the variance forecast starts at 0.409 and rises day by day towards the unconditional 0.58. The interval half-width grows with it; it ignores parameter uncertainty.

Limitations

• oe.arch accepts no weights.
• The common oe.predict and oe.margins adapters do not cover oe.arch; use its own result tables and family helpers.
• One series per fit, without gaps: missing="drop" removes rows only at the start or the end, and a missing value inside the series raises time_gaps.
• Because no positivity constraints are imposed, a lag the data do not support can make the likelihood unbounded; such fits, and power-ARCH fits whose power is 1 or less, are reported as nonconvergence with an explanation.
• Power ARCH is the symmetric form only; asymmetric power variants (such as Stata's aparch) and FIEGARCH are not implemented.
• The ARCH-in-mean term uses the current conditional variance only, with one of three fixed transformations (no archmlags()).
• EGARCH and power-ARCH variance forecasts beyond one step are plug-in forecasts of the transformed state, not conditional expectations, and all forecast standard errors ignore parameter uncertainty.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3), 307–327. doi:10.1016/0304-4076(86)90063-1
• Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987–1007. doi:10.2307/1912773
• 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. doi:10.1111/j.1540-6261.1993.tb05128.x
• Nelson, D. B. (1991). Conditional heteroskedasticity in asset returns: A new approach. Econometrica, 59(2), 347–370. doi:10.2307/2938260