diff --git a/arch/tests/univariate/test_censored_garch.py b/arch/tests/univariate/test_censored_garch.py new file mode 100644 index 0000000000..a9f8df7874 --- /dev/null +++ b/arch/tests/univariate/test_censored_garch.py @@ -0,0 +1,428 @@ +import numpy as np +from numpy.random import RandomState +from numpy.testing import assert_allclose, assert_equal +import pandas as pd +import pytest +from scipy.stats import norm + +from arch.univariate import ConstantMean +from arch.univariate.censored_garch import ( + CensoredGARCH, + _norm_pdf, + _norm_pdf_python, + _norm_sf, + _norm_sf_python, + censored_garch_recursion, + censored_garch_recursion_python, +) +from arch.univariate.volatility import GARCH + + +@pytest.fixture(scope="module") +def setup(): + rng = RandomState(1234) + t = 1000 + resids = rng.standard_normal(t) + threshold = 1.5 + censored = np.abs(resids) >= threshold + # Observed residuals are truncated at the threshold, as they would be + # for an actual censored return series. + observed = np.where(censored, np.sign(resids) * threshold, resids) + return { + "rng": rng, + "t": t, + "resids": observed, + "censored": censored, + "threshold": threshold, + "sigma2": np.zeros_like(observed), + } + + +def test_bounds_and_names(setup): + vol = CensoredGARCH(censored=setup["censored"], threshold=setup["threshold"]) + resids = setup["resids"] + + bounds = vol.bounds(resids) + v = np.mean(resids**2.0) + assert_allclose(bounds[0], (1.0e-8 * v, 10.0 * v)) + assert_equal(bounds[1], (0.0, 1.0)) + assert_equal(bounds[2], (0.0, 1.0)) + + assert vol.parameter_names() == ["omega", "alpha[1]", "beta[1]"] + assert vol.num_params == 3 + + +def test_constraints(setup): + vol = CensoredGARCH(censored=setup["censored"], threshold=setup["threshold"]) + a, b = vol.constraints() + # omega >= 0, alpha >= 0, beta >= 0, alpha + beta <= 1 + a_target = np.vstack((np.eye(3), np.array([[0, -1.0, -1.0]]))) + b_target = np.array([0.0, 0.0, 0.0, -1.0]) + assert_allclose(a, a_target) + assert_allclose(b, b_target) + + +def test_invalid_orders(): + censored = np.zeros(10, dtype=bool) + with pytest.raises(ValueError, match="One of p or q must be strictly positive"): + CensoredGARCH(censored=censored, threshold=1.0, p=0, q=0) + + +def test_invalid_threshold(): + censored = np.array([False, True, False]) + with pytest.raises(ValueError, match="threshold must be strictly positive"): + CensoredGARCH(censored=censored, threshold=np.array([1.0, -1.0, 1.0])) + + +def test_mismatched_threshold_length(): + censored = np.array([False, True, False]) + with pytest.raises(ValueError, match="same length"): + CensoredGARCH(censored=censored, threshold=np.array([1.0, 1.0])) + + +def test_no_censoring_matches_plain_garch(setup): + """With censored all-False, the recursion must reduce exactly to a plain + GARCH(1, 1) recursion, since the correction term is never applied.""" + + resids = setup["resids"] + t = setup["t"] + censored = np.zeros(t, dtype=bool) + + vol = CensoredGARCH(censored=censored, threshold=1.0, p=1, q=1) + vol._start, vol._stop = 0, t + backcast = vol.backcast(resids) + var_bounds = vol.variance_bounds(resids) + parameters = np.array([0.1, 0.1, 0.8]) + + sigma2_censored = np.zeros(t) + vol.compute_variance(parameters, resids, sigma2_censored, backcast, var_bounds) + + garch = GARCH(p=1, q=1) + sigma2_plain = np.zeros(t) + garch.compute_variance(parameters, resids, sigma2_plain, backcast, var_bounds) + + assert_allclose(sigma2_censored, sigma2_plain) + + +def test_censoring_increases_variance(setup): + """A censored observation's effective squared residual must be at least + as large as the naive (truncated) squared residual -- the whole point of + the correction -- and strictly larger whenever the threshold truly binds + below the model's own conditional standard deviation.""" + resids = setup["resids"] + censored = setup["censored"] + threshold = setup["threshold"] + t = setup["t"] + + vol_censored = CensoredGARCH(censored=censored, threshold=threshold, p=1, q=1) + vol_censored._start, vol_censored._stop = 0, t + vol_naive = CensoredGARCH( + censored=np.zeros(t, dtype=bool), threshold=threshold, p=1, q=1 + ) + vol_naive._start, vol_naive._stop = 0, t + + backcast = vol_censored.backcast(resids) + var_bounds = vol_censored.variance_bounds(resids) + parameters = np.array([0.05, 0.1, 0.85]) + + sigma2_censored = np.zeros(t) + vol_censored.compute_variance( + parameters, resids, sigma2_censored, backcast, var_bounds + ) + sigma2_naive = np.zeros(t) + vol_naive.compute_variance(parameters, resids, sigma2_naive, backcast, var_bounds) + + # sigma2 at t depends on sigma2/eff2 at t-1, so compare the *next* + # period's variance following a censored observation. + censored_idx = np.where(censored[:-1])[0] + assert len(censored_idx) > 0 + assert np.all(sigma2_censored[censored_idx + 1] >= sigma2_naive[censored_idx + 1]) + assert np.any(sigma2_censored[censored_idx + 1] > sigma2_naive[censored_idx + 1]) + + +def test_recursion_deep_tail_underflow_branch(): + """When the threshold is many standard deviations out, 1 - Phi(a) + underflows to 0 in floating point before a itself does. The recursion + should fall back to the asymptotic correction 1 + a^2 rather than + dividing by (effectively) zero.""" + + a = 40.0 # norm.sf(40) underflows to exactly 0.0 in float64 + assert norm.sf(a) == 0.0 + + threshold = a + # omega=0, alpha=1, backcast=1.0 => sigma2[0] = 0 + 1*backcast = 1.0 + # exactly, so a = threshold / sigma_t = threshold / 1.0 = a as intended. + parameters = np.array([0.0, 1.0, 0.0]) + resids = np.array([threshold, 0.0]) + censored = np.array([True, False]) + thresh_arr = np.array([threshold, 0.0]) + var_bounds = np.array([[1e-12, 1e12], [1e-12, 1e12]]) + + # Exercise both the (possibly jitted) recursion and the pure-Python source. + for recursion in (censored_garch_recursion, censored_garch_recursion_python): + sigma2 = np.zeros(2) + recursion( + parameters, resids, censored, thresh_arr, sigma2, 1, 1, 2, 1.0, var_bounds + ) + assert_allclose(sigma2[0], 1.0) + # sigma2[1] = alpha * eff2[0] = 1.0 * (1 + a^2) * sigma2[0] = 1 + a^2 + assert_allclose(sigma2[1], 1.0 + a * a, rtol=1e-10) + + +def test_recursion_zero_threshold_branch(): + """A censored observation with a zero threshold (a degenerate edge + case not reachable through the public API, which requires threshold > + 0 wherever censored is True) should fall back to the uncorrected + (corr=1.0) branch rather than raising or dividing by zero.""" + + # omega=0, alpha=1, backcast=1.0 => sigma2[0] = 1.0 exactly (isolates the + # step cleanly, same trick as the deep-tail test above). + parameters = np.array([0.0, 1.0, 0.0]) + resids = np.array([0.0, 0.0]) + censored = np.array([True, False]) + threshold = np.array([0.0, 0.0]) + var_bounds = np.array([[1e-12, 1e12], [1e-12, 1e12]]) + + # Exercise both the (possibly jitted) recursion and the pure-Python source. + for recursion in (censored_garch_recursion, censored_garch_recursion_python): + sigma2 = np.zeros(2) + recursion( + parameters, resids, censored, threshold, sigma2, 1, 1, 2, 1.0, var_bounds + ) + assert_allclose(sigma2[0], 1.0) + # corr = 1.0 (a <= 0 branch), so eff2[0] = sigma2[0] = 1.0 and + # sigma2[1] = alpha * eff2[0] = 1.0 + assert_allclose(sigma2[1], 1.0) + + +def test_recursion_zero_variance_branch(): + """When the model's own conditional variance is driven to exactly the + (zero) lower bound at a censored step, sigma_t is 0 and a would + otherwise be a division by zero; the recursion should fall back to + a = 0 (corr = 1.0) instead.""" + + parameters = np.array([0.0, 0.0, 0.0]) # omega=0 -> sigma2[0] forced to 0 + resids = np.array([1.5]) + censored = np.array([True]) + threshold = np.array([1.5]) + sigma2 = np.zeros(1) + var_bounds = np.array([[0.0, 1e12]]) # lower bound of exactly 0 + + censored_garch_recursion( + parameters, resids, censored, threshold, sigma2, 1, 1, 1, 0.0, var_bounds + ) + assert sigma2[0] == 0.0 + + +def test_norm_helpers_match_scipy(): + """The hand-rolled, numba-jittable _norm_pdf/_norm_sf (built from + math.exp/math.erfc rather than scipy.stats.norm, so they can be + inlined into the jitted recursion) must agree with scipy across a + wide range of inputs, including the deep tail where naive 1 - Phi(x) + would underflow long before erfc-based _norm_sf does.""" + + xs = np.array([0.0, 0.1, 0.5, 1.0, 1.96, 3.0, 5.0, 8.0, 15.0, 30.0]) + for x in xs: + # Check both the (possibly jitted) helpers and their pure-Python source. + for pdf in (_norm_pdf, _norm_pdf_python): + assert_allclose(pdf(x), norm.pdf(x), rtol=1e-12, atol=1e-300) + # scipy's own norm.sf underflows well before erfc-based _norm_sf + # does, so only compare where scipy itself hasn't hit zero. + scipy_sf = norm.sf(x) + if scipy_sf > 0: + for sf in (_norm_sf, _norm_sf_python): + assert_allclose(sf(x), scipy_sf, rtol=1e-8) + + +def test_jitted_recursion_matches_pure_python(): + """The numba-jitted censored_garch_recursion (used by compute_variance) + must produce bit-identical output to the pure-Python source it's + compiled from, censored_garch_recursion_python, across a realistic + GARCH(1,1)-then-censored path.""" + + rng = RandomState(7) + nobs = 2000 + threshold = 1.8 + resids = rng.standard_normal(nobs) + censored = np.abs(resids) >= threshold + resids = np.where(censored, np.sign(resids) * threshold, resids) + thresh_arr = threshold * np.ones(nobs) * censored + parameters = np.array([0.05, 0.1, 0.85]) + var_bounds = np.tile([1e-12, 1e12], (nobs, 1)) + + sigma2_jit = np.zeros(nobs) + sigma2_py = np.zeros(nobs) + censored_garch_recursion( + parameters, + resids, + censored, + thresh_arr, + sigma2_jit, + 1, + 1, + nobs, + 1.0, + var_bounds, + ) + censored_garch_recursion_python( + parameters, resids, censored, thresh_arr, sigma2_py, 1, 1, nobs, 1.0, var_bounds + ) + assert_allclose(sigma2_jit, sigma2_py, rtol=0, atol=0) + + +def test_correction_matches_truncated_normal_identity(): + """Directly check the recursion's censoring correction against the + truncated-normal second-moment identity + E[Z^2 | |Z| >= a] = 1 + a * phi(a) / (1 - Phi(a)) + for a single censored step with a known sigma2 carried in from omega.""" + + omega = 1.0 + threshold = 1.2 + parameters = np.array([omega, 0.0, 0.0]) # p=1, q=1, both zeroed out + resids = np.array([threshold]) + censored = np.array([True]) + thresh_arr = np.array([threshold]) + sigma2 = np.zeros(1) + var_bounds = np.array([[1e-12, 1e12]]) + + censored_garch_recursion( + parameters, resids, censored, thresh_arr, sigma2, 1, 1, 1, 1.0, var_bounds + ) + # sigma2[0] = omega = 1.0 regardless of censoring (alpha=beta=0 here); + # the correction only shows up in the *next* period's effective + # variance input, which we can't observe directly from sigma2 alone in + # this minimal example, so instead verify the corr factor formula used + # inside the recursion analytically. + a = threshold / np.sqrt(sigma2[0]) + expected_corr = 1.0 + a * norm.pdf(a) / norm.sf(a) + + # Re-run with a second observation to pull the correction into sigma2. + parameters2 = np.array([0.0, 1.0, 0.0]) # sigma2[1] = alpha * eff2[0] + resids2 = np.array([threshold, 0.0]) + censored2 = np.array([True, False]) + thresh2 = np.array([threshold, 0.0]) + sigma2_2 = np.zeros(2) + var_bounds2 = np.array([[1e-12, 1e12], [1e-12, 1e12]]) + censored_garch_recursion( + parameters2, resids2, censored2, thresh2, sigma2_2, 1, 1, 2, 1.0, var_bounds2 + ) + assert_allclose(sigma2_2[1], expected_corr, rtol=1e-10) + + +def test_fit_and_forecast_end_to_end(): + """Smoke test the full estimation pipeline: ConstantMean + CensoredGARCH + should fit without error and produce a finite one-step forecast.""" + + rng = RandomState(42) + nobs = 500 + threshold = 1.8 + resids = rng.standard_normal(nobs) + censored = np.abs(resids) >= threshold + observed = np.where(censored, np.sign(resids) * threshold, resids) + y = pd.Series(observed) + + model = ConstantMean(y) + model.volatility = CensoredGARCH(censored=censored, threshold=threshold, p=1, q=1) + res = model.fit(disp="off") + + assert np.all(np.isfinite(res.params.values)) + assert res.params["alpha[1]"] >= 0 + assert res.params["beta[1]"] >= 0 + + forecast = res.forecast(horizon=1, reindex=False) + assert np.isfinite(forecast.variance.values[-1, 0]) + assert forecast.variance.values[-1, 0] > 0 + + +def test_simulate_shapes(): + censored = np.zeros(10, dtype=bool) + vol = CensoredGARCH(censored=censored, threshold=1.0, p=1, q=1) + rng = RandomState(0) + parameters = np.array([0.05, 0.1, 0.85]) + data, sigma2 = vol.simulate(parameters, 250, rng.standard_normal, burn=50) + assert data.shape == (250,) + assert sigma2.shape == (250,) + assert np.all(sigma2 > 0) + + +def test_unsupported_forecast_options(): + censored = np.zeros(10, dtype=bool) + vol = CensoredGARCH(censored=censored, threshold=1.0, p=1, q=1) + with pytest.raises(NotImplementedError, match="horizon=1"): + vol._check_forecasting_method("analytic", horizon=2) + with pytest.raises(NotImplementedError, match="analytic"): + vol._check_forecasting_method("simulation", horizon=1) + with pytest.raises(NotImplementedError, match="Simulation-based forecasts"): + vol._simulation_forecast( + np.array([0.05, 0.1, 0.85]), + np.zeros(10), + 1.0, + np.tile([1e-8, 1e8], (10, 1)), + 0, + 1, + 10, + RandomState(0).standard_normal, + ) + + +def test_compute_variance_length_mismatch(): + """A censored/threshold series that doesn't line up with the data being + fit (neither nobs nor nobs - 1) must raise rather than silently misalign.""" + vol = CensoredGARCH(censored=np.zeros(10, dtype=bool), threshold=1.0) + vol._start, vol._stop = 0, 10 + resids = np.zeros(5) + var_bounds = np.tile([1e-8, 1e8], (5, 1)) + with pytest.raises(ValueError, match="does not match the data"): + vol.compute_variance( + np.array([0.05, 0.1, 0.85]), resids, np.zeros(5), 1.0, var_bounds + ) + + +@pytest.mark.parametrize(("p", "q"), [(1, 0), (0, 1), (2, 2)]) +def test_starting_values_orders(setup, p, q): + """starting_values must handle models with no ARCH term (p=0) or no + GARCH term (q=0), as well as higher orders.""" + resids = setup["resids"] + t = setup["t"] + vol = CensoredGARCH( + censored=setup["censored"], threshold=setup["threshold"], p=p, q=q + ) + vol._start, vol._stop = 0, t + sv = vol.starting_values(resids) + assert sv.shape == (1 + p + q,) + assert np.all(np.isfinite(sv)) + assert sv[0] > 0 + assert np.all(sv[1:] >= 0) + + +def test_simulate_initial_value(): + """A user-supplied initial_value is used as-is for the first variance.""" + vol = CensoredGARCH(censored=np.zeros(10, dtype=bool), threshold=1.0) + rng = RandomState(0) + parameters = np.array([0.05, 0.1, 0.85]) + data, sigma2 = vol.simulate( + parameters, 20, rng.standard_normal, burn=0, initial_value=0.5 + ) + assert data.shape == (20,) + assert_allclose(sigma2[0], 0.5) + + +def test_simulate_default_initial_value_stationary(): + """With persistence < 1 the default start is the unconditional variance.""" + vol = CensoredGARCH(censored=np.zeros(10, dtype=bool), threshold=1.0) + rng = RandomState(0) + parameters = np.array([0.05, 0.1, 0.85]) + _, sigma2 = vol.simulate(parameters, 20, rng.standard_normal, burn=0) + assert_allclose(sigma2[0], 0.05 / (1.0 - 0.95)) + + +def test_simulate_default_initial_value_unit_root(): + """With persistence >= 1 there is no unconditional variance, so the + default start falls back to omega.""" + vol = CensoredGARCH(censored=np.zeros(10, dtype=bool), threshold=1.0) + rng = RandomState(0) + parameters = np.array([0.05, 0.5, 0.5]) + data, sigma2 = vol.simulate(parameters, 20, rng.standard_normal, burn=0) + assert_allclose(sigma2[0], 0.05) + assert np.all(np.isfinite(data)) + assert np.all(sigma2 > 0) diff --git a/arch/univariate/__init__.py b/arch/univariate/__init__.py index 836f7fd0ac..230de2f342 100644 --- a/arch/univariate/__init__.py +++ b/arch/univariate/__init__.py @@ -1,6 +1,7 @@ import types from arch.univariate import recursions_python +from arch.univariate.censored_garch import CensoredGARCH from arch.univariate.diagnostics import excess_kurtosis, hill_estimator, var_ratio from arch.univariate.distribution import ( Distribution, @@ -44,6 +45,7 @@ "APARCH", "ARCH", "ARX", + "CensoredGARCH", "EGARCH", "FIAPARCH", "FIGARCH", diff --git a/arch/univariate/censored_garch.py b/arch/univariate/censored_garch.py new file mode 100644 index 0000000000..df4494095e --- /dev/null +++ b/arch/univariate/censored_garch.py @@ -0,0 +1,486 @@ +""" +Censored/latent-variance GARCH(p, q). + +Implements the variance-recursion correction described in Morgan & Trevor +(1999) and used in subsequent studies of censored return series (e.g. prices +subject to an exchange circuit-breaker or price-limit band): when the return +at time t is known to be censored (the true, latent shock would have exceeded +a known threshold in absolute value, but only the threshold itself is +observed), using the raw observed (truncated) squared residual in the GARCH +recursion understates the true conditional variance, because it discards the +information that the *true* shock was larger than what was recorded. + +This class corrects the recursion by replacing the squared residual at a +censored observation with its conditional expectation given censoring, under +the model's own (currently Gaussian) innovation distribution: + +.. math:: + + E[\\varepsilon_t^2 \\mid |\\varepsilon_t| \\ge c_t, \\mathcal{F}_{t-1}] + = \\sigma_t^2 \\left(1 + a_t \\cdot \\frac{\\phi(a_t)}{1 - \\Phi(a_t)}\\right), + \\quad a_t = c_t / \\sigma_t + +where :math:`c_t` is the (known) censoring threshold in effect at time t, +:math:`\\phi` and :math:`\\Phi` are the standard normal pdf/cdf, and +:math:`\\sigma_t^2` is the model's own conditional variance at t (available +within the recursion since it is computed sequentially forward). This is the +standard truncated-normal second-moment identity applied to a two-sided +(symmetric) censoring rule. + +Two importance caveats, stated in Notes on the class itself: + +* The correction uses the *model's own* :math:`\\sigma_t` for uncensored + and censored steps alike; it therefore self-consistently propagates the + extra variance from censored periods into subsequent, uncensored ones + through the ordinary GARCH recursion. +* The correction is derived for Gaussian innovations. The recursion is + still usable, as an approximation, with a non-Gaussian ``Distribution`` + (e.g. Student's t) selected for the *overall* model likelihood, but the + correction factor itself is not re-derived for that case. + +See the accompanying paper referenced in the class docstring for the +underlying empirical motivation and validation. +""" + +from arch.compat.numba import jit + +from collections.abc import Sequence +import math +from typing import cast + +import numpy as np + +from arch._typing import ( + ArrayLike1D, + Float64Array, + Float64Array1D, + Float64Array2D, + ForecastingMethod, + RNGType, +) +from arch.univariate.recursions_python import bounds_check +from arch.univariate.volatility import VarianceForecast, VolatilityProcess +from arch.utility.array import AbstractDocStringInheritor, ensure1d, to_array_1d + +__all__ = [ + "CensoredGARCH", + "censored_garch_recursion", + "censored_garch_recursion_python", +] + +_INV_SQRT_2PI = 1.0 / math.sqrt(2.0 * math.pi) +_SQRT2 = math.sqrt(2.0) + + +def _norm_pdf_python(x: float) -> float: + """Standard normal pdf, written with only numba-jittable primitives + (no scipy dependency) so this can be inlined into a jitted recursion.""" + return _INV_SQRT_2PI * math.exp(-0.5 * x * x) + + +def _norm_sf_python(x: float) -> float: + """Standard normal survival function 1 - Phi(x), via erfc rather than + 1 - erf, for the same numerical-stability reason scipy uses erfc + internally: it avoids catastrophic cancellation for large x, so this + underflows to exactly 0.0 only much further into the tail than a naive + 1 - Phi(x) would.""" + return 0.5 * math.erfc(x / _SQRT2) + + +# Jitted eagerly (rather than left to numba's lazy inlining) so these are +# ordinary numba Dispatcher objects the recursion below can call directly +# in nopython mode -- the same pattern used for `bounds_check` throughout +# arch.univariate.recursions_python. +_norm_pdf = jit(_norm_pdf_python, nopython=True, inline="always") +_norm_sf = jit(_norm_sf_python, nopython=True, inline="always") + + +def censored_garch_recursion_python( + parameters: Float64Array1D, + resids: Float64Array1D, + censored: Float64Array1D, + threshold: Float64Array1D, + sigma2: Float64Array1D, + p: int, + q: int, + nobs: int, + backcast: float, + var_bounds: Float64Array2D, +) -> Float64Array1D: + """ + Reference recursion for :class:`CensoredGARCH`. + + This is the pure-Python source; ``censored_garch_recursion`` (no + ``_python`` suffix) is this same function JIT-compiled with numba when + numba is installed, following the pattern used throughout + :mod:`arch.univariate.recursions_python`. When numba is unavailable, + ``arch.compat.numba.jit`` falls back to calling this function directly + (with a one-time performance warning), so correctness never depends on + numba being installed -- only speed does. + + Parameters + ---------- + parameters : ndarray + Model parameters, ordered [omega, alpha_1..alpha_p, beta_1..beta_q]. + resids : ndarray + Mean-equation residuals. For a censored observation, this is the + *observed* (truncated) residual, i.e. it equals +/- the threshold + in effect at that time, not the (unobserved) true shock. + censored : ndarray + Same length as resids. Nonzero (truthy) where the observation at + that index is censored. + threshold : ndarray + Same length as resids. The censoring threshold |c_t| in effect at + each time t (only used where ``censored`` is truthy). + sigma2 : ndarray + Output array for the conditional variance. + p, q : int + GARCH orders (no asymmetric/leverage term in this first version). + nobs : int + Number of observations. + backcast : float + Initial value used before any real observations are available. + var_bounds : ndarray + Per-observation [lower, upper] bounds enforced on sigma2. + + Returns + ------- + sigma2 : ndarray + The conditional variance, written in place and also returned. + """ + eff2 = np.empty(nobs, dtype=float) + + for t in range(nobs): + loc = 0 + sigma2[t] = parameters[loc] + loc += 1 + for j in range(p): + idx = t - 1 - j + if idx < 0: + sigma2[t] += parameters[loc] * backcast + else: + sigma2[t] += parameters[loc] * eff2[idx] + loc += 1 + for j in range(q): + idx = t - 1 - j + if idx < 0: + sigma2[t] += parameters[loc] * backcast + else: + sigma2[t] += parameters[loc] * sigma2[idx] + loc += 1 + + sigma2[t] = bounds_check(sigma2[t], var_bounds[t]) + + if not censored[t]: + eff2[t] = resids[t] * resids[t] + else: + sigma_t = np.sqrt(sigma2[t]) + a = threshold[t] / sigma_t if sigma_t > 0 else 0.0 + if a <= 0: + corr = 1.0 + else: + one_minus_phi = _norm_sf(a) + if one_minus_phi < 1e-12: + # Far into the tail: 1 - Phi(a) underflows before a does. + # Asymptotically phi(a)/(1-Phi(a)) ~ a, so the correction + # factor 1 + a*hazard(a) -> 1 + a^2. + corr = 1.0 + a * a + else: + corr = 1.0 + a * _norm_pdf(a) / one_minus_phi + eff2[t] = corr * sigma2[t] + + return sigma2 + + +# JIT-compiled with numba when available (falls back to the pure-Python +# function above, with a one-time performance warning, when it isn't). +censored_garch_recursion = jit(censored_garch_recursion_python, nopython=True) + + +class CensoredGARCH(VolatilityProcess, metaclass=AbstractDocStringInheritor): + r""" + Latent-variance GARCH(p, q) for censored return series. + + Models return series subject to a known, time-varying censoring rule + (e.g. an exchange price-limit band). + + Parameters + ---------- + censored : {ndarray, Series} + Boolean (or 0/1) array, same length as the data used to fit the + model, indicating which observations are censored (the recorded + return equals the threshold in effect that period, rather than the + true underlying shock). + threshold : {ndarray, Series, float} + The censoring threshold |c_t| in effect at each time t, in the same + units as the model residuals. May be passed as a scalar if the + threshold is constant over the sample. + p : int, optional + Order of the symmetric innovation. Default is 1. + q : int, optional + Order of the lagged conditional variance. Default is 1. + + Notes + ----- + In this class of processes, the variance dynamics are + + .. math:: + + \sigma_t^2 = \omega + \sum_{i=1}^p \alpha_i \tilde\varepsilon_{t-i}^2 + + \sum_{k=1}^q \beta_k \sigma_{t-k}^2 + + where :math:`\tilde\varepsilon_s^2 = \varepsilon_s^2` when observation + :math:`s` is not censored, and + + .. math:: + + \tilde\varepsilon_s^2 = \sigma_s^2 + \left(1 + a_s \frac{\phi(a_s)}{1-\Phi(a_s)}\right), \quad + a_s = c_s / \sigma_s + + when it is (the conditional expectation of the squared shock given that + it exceeded the threshold, under Gaussian innovations). This uses the + truncated-normal second-moment identity + :math:`\int_a^\infty x^2\phi(x)\,dx = a\phi(a) + (1-\Phi(a))`. + + Only symmetric (no leverage/asymmetric term) GARCH(p, q) is supported in + this first version; ``power`` is fixed at 2.0. The recursion + (:func:`censored_garch_recursion`) is JIT-compiled with numba when numba + is installed -- about 100x faster than the pure-Python source it is + compiled from (:func:`censored_garch_recursion_python`) in informal + benchmarking -- following the same pattern used throughout + :mod:`arch.univariate.recursions_python`; correctness does not depend on + numba being installed, only speed does. There is no separate Cython + backend. + + Analytic multi-step forecasting and simulation-based forecasting are not + yet implemented (only one-step-ahead analytic forecasts, and forward + simulation assuming no future censoring, are supported). Extending + analytic forecasting to horizon > 1 is not a simple drop-in of the + approach used by :class:`GARCH`: doing it correctly requires threading + the corrected effective squared residual (``eff2`` inside the + recursion, not the raw truncated residual) through the lag lookback for + any censored observation within the forecast origin's lag window, which + the recursion does not currently expose. + + Examples + -------- + >>> from arch.univariate import CensoredGARCH + >>> cens = (returns.abs() >= 0.0995) # e.g. a +/-10% price-limit band + >>> vol = CensoredGARCH(censored=cens, threshold=0.0995) + """ + + def __init__( + self, + censored: ArrayLike1D, + threshold: ArrayLike1D | float, + p: int = 1, + q: int = 1, + ) -> None: + super().__init__() + if p < 0 or q < 0 or (p == 0 and q == 0): + raise ValueError("One of p or q must be strictly positive") + self.p: int = int(p) + self.q: int = int(q) + self._num_params = 1 + self.p + self.q + self._name = f"Censored-GARCH({self.p}, {self.q})" + self._updatable = False + + censored_arr = ensure1d(censored, "censored", True) + self._censored = np.asarray(censored_arr, dtype=bool) + + if np.isscalar(threshold): + self._threshold = np.full( + self._censored.shape[0], float(threshold), dtype=float + ) + else: + self._threshold = to_array_1d( + ensure1d(threshold, "threshold", True) + ).astype(float) + if self._threshold.shape[0] != self._censored.shape[0]: + raise ValueError( + "censored and threshold must be the same length when " + "threshold is not a scalar" + ) + if np.any(self._threshold[self._censored] <= 0): + raise ValueError( + "threshold must be strictly positive wherever censored is True" + ) + + def bounds(self, resids: ArrayLike1D) -> list[tuple[float, float]]: + v = float(np.mean(np.asarray(resids) ** 2)) + bounds = [(1e-8 * v, 10.0 * v)] + bounds.extend([(0.0, 1.0)] * self.p) + bounds.extend([(0.0, 1.0)] * self.q) + return bounds + + def constraints(self) -> tuple[Float64Array, Float64Array]: + k = self.p + self.q + a = np.zeros((k + 2, k + 1)) + for i in range(k + 1): + a[i, i] = 1.0 + a[k + 1, 1:] = -1.0 + b = np.zeros(k + 2) + b[k + 1] = -1.0 + return a, b + + def _slice(self, arr: Float64Array1D) -> Float64Array1D: + return arr[self._start : self._stop] + + def compute_variance( + self, + parameters: Float64Array1D, + resids: ArrayLike1D, + sigma2: Float64Array1D, + backcast: float | Float64Array1D, + var_bounds: Float64Array2D, + ) -> Float64Array1D: + _resids = to_array_1d(resids) + nobs = _resids.shape[0] + censored = self._slice(self._censored) + threshold = self._slice(self._threshold) + if censored.shape[0] == nobs - 1: + # `_one_step_forecast` (the base class helper used by + # `_analytic_forecast`) appends one extra placeholder residual + # (0.0) past the end of the fitted sample to obtain a one-step + # ahead variance forecast. That placeholder observation isn't + # part of the real censored/threshold series, and it is never + # itself fed back into the recursion (it is always the last + # point), so it is safe to pad with an "uncensored" placeholder. + censored = np.concatenate([censored, np.array([False])]) + threshold = np.concatenate([threshold, np.array([0.0])]) + elif censored.shape[0] != nobs: + raise ValueError( + "The length of the censored/threshold arrays passed at " + "construction does not match the data used to fit the " + "model. Re-check `start`/`stop` slicing." + ) + assert isinstance(backcast, float) + censored_garch_recursion( + parameters, + _resids, + censored, + threshold, + sigma2, + self.p, + self.q, + nobs, + backcast, + var_bounds, + ) + return sigma2 + + def starting_values(self, resids: ArrayLike1D) -> Float64Array1D: + p, q = self.p, self.q + target = float(np.mean(np.asarray(resids) ** 2)) + var_bounds = self.variance_bounds(resids) + backcast = self.backcast(resids) + + alphas = [0.05, 0.1, 0.2] + persistences = [0.7, 0.9, 0.98] + svs: list[Float64Array1D] = [] + llfs = [] + for alpha_total, persistence in [(a, b) for a in alphas for b in persistences]: + sv = np.zeros(1 + p + q, dtype=float) + sv[0] = (1.0 - persistence) * target + if p > 0: + sv[1 : 1 + p] = alpha_total / p + if q > 0: + sv[1 + p : 1 + p + q] = (persistence - alpha_total) / q + sv[1 + p : 1 + p + q] = np.clip(sv[1 + p : 1 + p + q], 1e-4, None) + svs.append(cast("Float64Array1D", sv)) + llfs.append(self._gaussian_loglikelihood(sv, resids, backcast, var_bounds)) + loc = int(np.argmax(llfs)) + return svs[loc] + + def parameter_names(self) -> list[str]: + names = ["omega"] + names.extend([f"alpha[{i + 1}]" for i in range(self.p)]) + names.extend([f"beta[{i + 1}]" for i in range(self.q)]) + return names + + def simulate( + self, + parameters: Sequence[int | float] | ArrayLike1D, + nobs: int, + rng: RNGType, + burn: int = 500, + initial_value: float | Float64Array | None = None, + ) -> tuple[Float64Array, Float64Array]: + # Simulates from the model assuming NO censoring in the simulated path + # (future censoring status is not known ex ante). This is the + # *unconditional* process the model implies: exact only when the true + # return process is never censored, and an approximation whenever a + # future censoring rule would actually bind. See the class Notes. + parameters = ensure1d(parameters, "parameters", False) + p, q = self.p, self.q + errors = rng(nobs + burn) + + if initial_value is None: + persistence = float(np.sum(parameters[1:])) + if (1.0 - persistence) > 0: + initial_value = parameters[0] / (1.0 - persistence) + else: + initial_value = parameters[0] + + sigma2 = np.zeros(nobs + burn) + data = np.zeros(nobs + burn) + max_lag = max(p, q, 1) + sigma2[:max_lag] = initial_value + data[:max_lag] = np.sqrt(sigma2[:max_lag]) * errors[:max_lag] + + for t in range(max_lag, nobs + burn): + loc = 0 + sigma2[t] = parameters[loc] + loc += 1 + for j in range(p): + sigma2[t] += parameters[loc] * data[t - 1 - j] ** 2 + loc += 1 + for j in range(q): + sigma2[t] += parameters[loc] * sigma2[t - 1 - j] + loc += 1 + data[t] = errors[t] * np.sqrt(sigma2[t]) + + return data[burn:], sigma2[burn:] + + def _check_forecasting_method( + self, method: ForecastingMethod, horizon: int + ) -> None: + if horizon > 1: + raise NotImplementedError( + "Multi-step forecasts are not yet implemented for " + "CensoredGARCH; only horizon=1 is supported." + ) + if method != "analytic": + raise NotImplementedError( + "Only method='analytic' is currently supported for CensoredGARCH." + ) + + def _analytic_forecast( + self, + parameters: Float64Array1D, + resids: ArrayLike1D, + backcast: float | Float64Array1D, + var_bounds: Float64Array2D, + start: int, + horizon: int, + ) -> VarianceForecast: + _, forecasts = self._one_step_forecast( + parameters, to_array_1d(resids), backcast, var_bounds, horizon, start + ) + return VarianceForecast(forecasts) + + def _simulation_forecast( + self, + parameters: Float64Array1D, + resids: ArrayLike1D, + backcast: float | Float64Array1D, + var_bounds: Float64Array2D, + start: int, + horizon: int, + simulations: int, + rng: RNGType, + ) -> VarianceForecast: + raise NotImplementedError( + "Simulation-based forecasts are not yet implemented for CensoredGARCH." + ) diff --git a/doc/source/changes/8.0.rst b/doc/source/changes/8.0.rst index aa3c86fca3..b7f5750e35 100644 --- a/doc/source/changes/8.0.rst +++ b/doc/source/changes/8.0.rst @@ -2,6 +2,15 @@ Version 8 ========= +Since Last Release +=================== + +- Added :class:`~arch.univariate.CensoredGARCH`, a GARCH(p, q) volatility + process for return series subject to a known, time-varying censoring rule + (e.g. an exchange price-limit band), using the truncated-normal + second-moment correction of Morgan and Trevor (1999). The recursion is + numba-JIT-compiled when numba is available. + Release 8.0 =========== - Moved to Meson build system from setuptools. diff --git a/doc/source/univariate/volatility.rst b/doc/source/univariate/volatility.rst index 97207dd8d1..df6954b100 100644 --- a/doc/source/univariate/volatility.rst +++ b/doc/source/univariate/volatility.rst @@ -22,6 +22,7 @@ volatility. ARCH APARCH FIAPARCH + CensoredGARCH Fractionally Integrated Asymmetric Power ARCH ---------------------------------------------