Volatility Persistence: Regimes vs Long Memory

The question

Every GARCH fit in the volatility arc reports the same near-unit-root persistence, $\alpha+\beta\approx 0.997$ — volatility shocks that almost never die. But why? Two rival stories explain it, and Diebold & Inoue (2001) proved they are observationally confounded — so this is a genuine contest. Long memory (FIGARCH): the persistence is real, a hyperbolically-decaying memory. Regime switching (MS-GARCH; Lamoureux–Lastrapes 1990): the persistence is spurious, an artifact of unmodelled shifts in the volatility level. This example fits both from scratch (each with an R package cross-check) and settles which one the S&P 500 supports.

FIGARCH — genuine long memory

FIGARCH (Baillie, Bollerslev & Mikkelsen 1996) offers the honest middle between short-memory GARCH ( $d=0$ ) and integrated IGARCH ( $d=1$ ): it applies the fractional-differencing operator $(1-L)^d$ to the ARCH filter, so a volatility shock decays hyperbolically ( $\sim j^{\,d-1}$ ) rather than geometrically. The estimate is $d=0.52$ — well inside $(0,1)$ , genuine long memory — and FIGARCH-t fits better than GARCH-t by AIC. Read as memory: GARCH has forgotten a shock within ~6 weeks, but FIGARCH still carries ~7% of it after 100 trading days. This is the same long memory the realized-volatility example found in intraday data — here recovered from daily returns alone.

\phi(L)\,(1-L)^d\,\varepsilon_t^2 = \omega + (1-\beta L)\,\nu_t, \qquad \nu_t = \varepsilon_t^2 - \sigma_t^2, \qquad d \in (0,1)

MS-GARCH — regime switching

MS-GARCH lets a latent 2-state Markov chain $S_t$ (calm / turbulent) switch the volatility dynamics. The trap: a naïve MS-GARCH makes $h_t$ depend on the entire regime path ( $2^t$ of them), so the likelihood is intractable. The Haas–Mittnik–Paolella (2004) specification gives each regime its own GARCH run on the full series, restoring a tractable likelihood via the Hamilton (1989) filter (smoothed regime probabilities from the Kim 1994 smoother). It needs Student-t innovations to find sustained regimes — a calm state (~0.66% daily vol) and a turbulent one (~1.45%) that track the dot-com bust and 2008 — rather than one-day tail spikes.

h_t^{(k)} = \omega_k + \alpha_k\, r_{t-1}^2 + \beta_k\, h_{t-1}^{(k)}, \qquad r_t \mid S_t{=}k \sim N\big(0, h_t^{(k)}\big)\ \text{or}\ t_\nu

The verdict

The verdict: both — but the persistence is long memory. The decisive fact: even though the regime model fits best — MS-GARCH-t posts the lowest AIC of the four — allowing regimes does not dissolve the persistence: the within-regime persistence stays at 0.99 (calm regime 0.991 vs the pooled 0.997). If the near-unit-root persistence were merely a switching artifact, accounting for the switches would have collapsed it — it didn't. And the realized-volatility evidence points the same way (a smooth hyperbolic autocorrelation, not discrete level shifts). So it is not either/or: the S&P 500 has both a slow regime structure (the vol level drifts with the business cycle) and genuine long memory (the persistence within each state). The near-unit-root persistence at the heart of the entire volatility arc is real, not a statistical illusion. Both from-scratch engines are confirmed by a package — FIGARCH via rugarch ( $d=0.52$ ), MS-GARCH via the MSGARCH package.

Notebooks

Python FIGARCH from scratch — fractional integration, d = 0.52, hyperbolic memory (Python) Open → R FIGARCH in R — rugarch fiGARCH cross-check (R) Open → Python MS-GARCH from scratch — Haas–Mittnik–Paolella regimes, Hamilton filter + Kim smoother (Python) Open → R MS-GARCH in R — MSGARCH package (per-regime GARCH, ML & Bayesian) (R) Open →

Downloads

msgarch.py Markov-switching GARCH (Haas–Mittnik–Paolella) — Hamilton filter + Kim smoother (NumPy / scipy) figarch.py FIGARCH — fractional-integration long-memory GARCH (NumPy / scipy) sp500ret.csv Daily S&P 500 returns, 1987–2009 spx_rv.csv S&P 500 realized variance — the long-memory cross-check in the MS-GARCH notebook

References

Baillie, R. T., Bollerslev, T. & Mikkelsen, H. O. (1996). Fractionally integrated generalized autoregressive conditional heteroskedasticity. Journal of Econometrics 74(1), 3–30. — the FIGARCH long-memory model
Diebold, F. X. & Inoue, A. (2001). Long memory and regime switching. Journal of Econometrics 105(1), 131–159. — the observational confounding of long memory and regimes
Lamoureux, C. G. & Lastrapes, W. D. (1990). Persistence in variance, structural change, and the GARCH model. Journal of Business & Economic Statistics 8(2), 225–234. — structural change as spurious GARCH persistence
Granger, C. W. J. & Hyung, N. (2004). Occasional structural breaks and long memory. Journal of Empirical Finance 11(3), 399–421. — breaks that mimic long memory in S&P 500 volatility
Haas, M., Mittnik, S. & Paolella, M. S. (2004). A new approach to Markov-switching GARCH models. Journal of Financial Econometrics 2(4), 493–530. — the per-regime MS-GARCH specification implemented here
Ardia, D., Bluteau, K., Boudt, K., Catania, L. & Trottier, D.-A. (2019). Markov-switching GARCH models in R: the MSGARCH package. Journal of Statistical Software 91(4), 1–38. — the R MSGARCH cross-check

MS-GARCH Module — Source Code

"""
Markov-switching GARCH (Haas, Mittnik & Paolella 2004 specification) from scratch, Normal or Student-t.

A latent 2-state Markov chain S_t (calm / turbulent) switches the volatility dynamics. In the HMP
specification each regime carries its OWN GARCH process, run on the full return series, which removes
the path-dependence that makes exact MS-GARCH likelihoods intractable:

    h_t^{(k)} = omega_k + alpha_k r_{t-1}^2 + beta_k h_{t-1}^{(k)},     k = 1,2
    r_t | S_t = k  ~  N(0, h_t^{(k)})  or  standardized t_nu with variance h_t^{(k)}
    S_t Markov with transition P[i,j] = P(S_t=j | S_{t-1}=i)

Likelihood via the Hamilton (1989) filter; smoothed regime probabilities via the Kim (1994) smoother.
theta = (o1,a1,b1, o2,a2,b2, p11, p22)  [+ nu if Student-t, one shared tail index across regimes].
"""
import numpy as np
from scipy.optimize import minimize
from scipy.stats import t as _t


def _paths(r2, th):
    o1, a1, b1, o2, a2, b2 = th[:6]
    T = len(r2); H = np.empty((T, 2))
    H[0, 0] = o1 / max(1 - a1 - b1, 1e-4); H[0, 1] = o2 / max(1 - a2 - b2, 1e-4)
    for t in range(1, T):
        H[t, 0] = o1 + a1 * r2[t - 1] + b1 * H[t - 1, 0]
        H[t, 1] = o2 + a2 * r2[t - 1] + b2 * H[t - 1, 1]
    return H


def _feasible(th, student):
    o1, a1, b1, o2, a2, b2, p11, p22 = th[:8]
    ok = (o1 > 0 and o2 > 0 and a1 >= 0 and b1 >= 0 and a2 >= 0 and b2 >= 0
          and a1 + b1 < 1 and a2 + b2 < 1 and 0 < p11 < 1 and 0 < p22 < 1)
    if student:
        ok = ok and th[8] > 2.05
    return ok


def _emissions(r, H, th, student):
    """Per-regime observation densities eta_t(k) = f(r_t | S_t=k), vectorised over t (shape T x 2)."""
    if student:
        nu = th[8]; sc = np.sqrt(H * (nu - 2.0) / nu)          # scale so Var(r_t|k)=h_t^{(k)}
        return _t.pdf(r[:, None], df=nu, scale=sc)
    return np.exp(-0.5 * (np.log(2 * np.pi * H) + r[:, None] ** 2 / H))


def filter(th, r, student=False, return_all=False):
    """Hamilton filter: log-likelihood (+ filtered probs, H paths, predicted probs if return_all)."""
    r2 = r ** 2; T = len(r); H = _paths(r2, th)
    p11, p22 = th[6], th[7]
    P = np.array([[p11, 1 - p11], [1 - p22, p22]])
    pi = np.array([1 - p22, 1 - p11]); pi = pi / pi.sum()        # stationary distribution
    eta = _emissions(r, H, th, student)
    xi = pi.copy(); ll = 0.0; filt = np.empty((T, 2)); pred = np.empty((T, 2))
    for t in range(T):
        pr = xi @ P                                             # P(S_t | data<t)
        num = pr * eta[t]; lik = num.sum()
        if lik <= 0 or not np.isfinite(lik):
            return (-np.inf, None, None, None) if return_all else -np.inf
        ll += np.log(lik); xi = num / lik
        filt[t] = xi; pred[t] = pr
    return (ll, filt, H, pred) if return_all else ll


def smooth(th, r, student=False):
    """Kim smoother: P(S_t = k | all data)."""
    _, filt, H, pred = filter(th, r, student, return_all=True)
    T = len(r); p11, p22 = th[6], th[7]
    P = np.array([[p11, 1 - p11], [1 - p22, p22]])
    sm = np.empty((T, 2)); sm[-1] = filt[-1]
    for t in range(T - 2, -1, -1):
        for j in range(2):
            sm[t, j] = filt[t, j] * np.sum(P[j, :] * sm[t + 1, :] / np.where(pred[t + 1] > 0, pred[t + 1], 1e-12))
        sm[t] /= sm[t].sum()
    return sm, H


def loglik(th, r, student=False):
    if not _feasible(th, student):
        return -np.inf
    return filter(th, r, student)


def mle(r, student=False):
    r = np.asarray(r, float); r = r - r.mean(); v = np.var(r)
    x0 = [0.02 * v, 0.05, 0.90, 0.30 * v, 0.15, 0.75, 0.98, 0.95] + ([8.0] if student else [])
    nll = lambda th: -loglik(th, r, student) if np.isfinite(loglik(th, r, student)) else 1e12
    res = minimize(nll, np.array(x0), method="Nelder-Mead", options=dict(xatol=1e-7, fatol=1e-7, maxiter=50000))
    th = res.x
    u1 = th[0] / (1 - th[1] - th[2]); u2 = th[3] / (1 - th[4] - th[5])
    if u1 > u2:                                                  # order so regime 2 is the HIGH-vol state
        tail = list(th[8:])
        th = np.array([th[3], th[4], th[5], th[0], th[1], th[2], th[7], th[6]] + tail)
    return th, r


def names(student=False):
    return ["omega1", "alpha1", "beta1", "omega2", "alpha2", "beta2", "p11", "p22"] + (["nu"] if student else [])