From 41b38464d7593950f6bdbd212984e00eb3ef3746 Mon Sep 17 00:00:00 2001 From: Muhammad Uzair Date: Mon, 28 Sep 2026 14:18:56 +0500 Subject: [PATCH 1/6] Add CensoredGARCH: latent-variance correction for censored return series Implements a GARCH(p, q) volatility process for return series subject to a known, time-varying censoring rule (e.g. an exchange price-limit band, such as PSX's daily circuit breaker). When a return is censored (the true shock exceeded a known threshold but only the threshold is observed), using the raw truncated squared residual in the ordinary GARCH recursion understates the true conditional variance. CensoredGARCH corrects this by replacing the squared residual at a censored observation with its conditional expectation given censoring, under the model's own conditional variance, using the truncated-normal second-moment identity: E[eps_t^2 | |eps_t| >= c_t] = sigma_t^2 * (1 + a_t * phi(a_t) / (1 - Phi(a_t))) a_t = c_t / sigma_t This is the correction described in Morgan & Trevor (1999) and applied to exchange price-limit censoring in follow-up work (Wei 2002, Hsieh & Yang 2009). A companion open-source package (mtgarch on PyPI) applies this same estimator to Pakistan Stock Exchange data; this PR ports the estimator into `arch` as a proper VolatilityProcess subclass so it can be used directly with arch_model's full estimation, inference, and forecasting machinery instead of a standalone implementation. - New `arch/univariate/censored_garch.py`: `CensoredGARCH` class and its pure-Python recursion (`censored_garch_recursion`), following the `GARCH`/`FixedVariance` classes as templates. Implements all required abstract methods (bounds, constraints, compute_variance, starting_values, parameter_names, simulate) plus one-step-ahead analytic forecasting. - `arch/univariate/__init__.py`: export `CensoredGARCH`. - `arch/tests/univariate/test_censored_garch.py`: unit tests covering parameter/constraint setup, input validation, the reduction to plain GARCH when nothing is censored, a direct check of the recursion against the truncated-normal identity, and an end-to-end ConstantMean + CensoredGARCH fit/forecast smoke test. - Docs: added to the volatility-processes autosummary and changelog. No accelerated (Cython/Numba) backend yet; only horizon-1 analytic forecasting is supported. Both are noted as known limitations in the class docstring and would be natural follow-ups. --- arch/tests/univariate/test_censored_garch.py | 228 ++++++++++ arch/univariate/__init__.py | 2 + arch/univariate/censored_garch.py | 436 +++++++++++++++++++ doc/source/changes/8.0.rst | 16 +- doc/source/univariate/volatility.rst | 1 + 5 files changed, 676 insertions(+), 7 deletions(-) create mode 100644 arch/tests/univariate/test_censored_garch.py create mode 100644 arch/univariate/censored_garch.py diff --git a/arch/tests/univariate/test_censored_garch.py b/arch/tests/univariate/test_censored_garch.py new file mode 100644 index 0000000000..6c49b9a471 --- /dev/null +++ b/arch/tests/univariate/test_censored_garch.py @@ -0,0 +1,228 @@ +import numpy as np +from numpy.random import RandomState +from numpy.testing import assert_allclose, assert_equal +import pytest +from scipy.stats import norm + +from arch.univariate import ConstantMean +from arch.univariate.censored_garch import CensoredGARCH + + +@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.""" + from arch.univariate.volatility import GARCH + + 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_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.""" + from arch.univariate.censored_garch import censored_garch_recursion + + 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.""" + import pandas as pd + + 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, + ) diff --git a/arch/univariate/__init__.py b/arch/univariate/__init__.py index 0f30e5980c..b5e16af65b 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.distribution import ( Distribution, GeneralizedError, @@ -41,6 +42,7 @@ "APARCH", "ARCH", "ARX", + "CensoredGARCH", "EGARCH", "FIGARCH", "GARCH", diff --git a/arch/univariate/censored_garch.py b/arch/univariate/censored_garch.py new file mode 100644 index 0000000000..d016d7030f --- /dev/null +++ b/arch/univariate/censored_garch.py @@ -0,0 +1,436 @@ +""" +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 collections.abc import Sequence +from typing import cast + +import numpy as np +from scipy.stats import norm + +from arch._typing import ( + ArrayLike1D, + Float64Array, + Float64Array1D, + Float64Array2D, + ForecastingMethod, + RNGType, +) +from arch.univariate.recursions_python import bounds_check_python as bounds_check +from arch.univariate.volatility import VarianceForecast, VolatilityProcess +from arch.utility.array import AbstractDocStringInheritor, ensure1d, to_array_1d + +__all__ = ["CensoredGARCH"] + + +def censored_garch_recursion( + parameters: Float64Array1D, + resids: Float64Array1D, + censored: Float64Array1D, + threshold: Float64Array1D, + sigma2: Float64Array1D, + p: int, + q: int, + nobs: int, + backcast: float, + var_bounds: Float64Array2D, +) -> Float64Array1D: + """ + Pure-Python reference recursion for :class:`CensoredGARCH`. + + Not currently accelerated (no Cython/Numba backend) -- see the class + docstring's "Notes" section. + + Parameters + ---------- + parameters : ndarray + [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 + + +class CensoredGARCH(VolatilityProcess, metaclass=AbstractDocStringInheritor): + r""" + Censored/latent-variance GARCH(p, q) for 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. + + 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) + + 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. No accelerated + (Cython/Numba) recursion is provided yet -- see + :func:`censored_garch_recursion`. + + 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). + """ + + 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]: + """ + Simulate from the model, assuming NO censoring in the simulated + path (future censoring status is not known ex ante -- this + simulates the *unconditional* process the model implies, which is + only exactly correct when the true return process is never + censored; it is an approximation whenever a future censoring rule + would actually bind). See 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: + sigma2, 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..813a898681 100644 --- a/doc/source/changes/8.0.rst +++ b/doc/source/changes/8.0.rst @@ -2,12 +2,14 @@ Version 8 ========= -Release 8.0 -=========== -- Moved to Meson build system from setuptools. -- Minimum required Python is now 3.10. -- Improved compatability with future NumPy and pandas, including pandas 3. -- Many small documentation fixes. -- Improved typing accuracy across the codebase. +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). +Release 8.0 +=========== +- Moved to Meson build system from setuptools. \ No newline at end of file diff --git a/doc/source/univariate/volatility.rst b/doc/source/univariate/volatility.rst index b3bda68573..c18f153a45 100644 --- a/doc/source/univariate/volatility.rst +++ b/doc/source/univariate/volatility.rst @@ -21,6 +21,7 @@ volatility. MIDASHyperbolic ARCH APARCH + CensoredGARCH Parameterless Variance Processes -------------------------------- From 443af255f874eb879d0a47e6be727c2c2c6de7c4 Mon Sep 17 00:00:00 2001 From: Muhammad Uzair Date: Mon, 28 Sep 2026 14:54:04 +0500 Subject: [PATCH 2/6] Add tests for censored_garch_recursion edge-case branches flagged by Codecov --- arch/tests/univariate/test_censored_garch.py | 74 ++++++++++++++++++++ 1 file changed, 74 insertions(+) diff --git a/arch/tests/univariate/test_censored_garch.py b/arch/tests/univariate/test_censored_garch.py index 6c49b9a471..e36617b906 100644 --- a/arch/tests/univariate/test_censored_garch.py +++ b/arch/tests/univariate/test_censored_garch.py @@ -131,6 +131,80 @@ def test_censoring_increases_variance(setup): 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.""" + from arch.univariate.censored_garch import censored_garch_recursion + from scipy.stats import norm + + 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]) + sigma2 = np.zeros(2) + var_bounds = np.array([[1e-12, 1e12], [1e-12, 1e12]]) + + censored_garch_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.""" + from arch.univariate.censored_garch import censored_garch_recursion + + # 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]) + sigma2 = np.zeros(2) + var_bounds = np.array([[1e-12, 1e12], [1e-12, 1e12]]) + + censored_garch_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.""" + from arch.univariate.censored_garch import censored_garch_recursion + + 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_correction_matches_truncated_normal_identity(): """Directly check the recursion's censoring correction against the truncated-normal second-moment identity From 19ed5a29eacf3f46ab39bd334d4d35091c56445d Mon Sep 17 00:00:00 2001 From: Muhammad Uzair Date: Mon, 28 Sep 2026 16:39:32 +0500 Subject: [PATCH 3/6] Accelerate CensoredGARCH recursion with numba JIT (~100x) --- arch/tests/univariate/test_censored_garch.py | 49 +++++++++++++ arch/univariate/censored_garch.py | 74 ++++++++++++++++---- doc/source/changes/8.0.rst | 11 ++- 3 files changed, 119 insertions(+), 15 deletions(-) diff --git a/arch/tests/univariate/test_censored_garch.py b/arch/tests/univariate/test_censored_garch.py index e36617b906..7790861a9e 100644 --- a/arch/tests/univariate/test_censored_garch.py +++ b/arch/tests/univariate/test_censored_garch.py @@ -205,6 +205,55 @@ def test_recursion_zero_variance_branch(): 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.""" + from arch.univariate.censored_garch import _norm_pdf, _norm_sf + + 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: + assert_allclose(_norm_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: + assert_allclose(_norm_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.""" + from arch.univariate.censored_garch import ( + censored_garch_recursion, + censored_garch_recursion_python, + ) + + 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 diff --git a/arch/univariate/censored_garch.py b/arch/univariate/censored_garch.py index d016d7030f..881163ab5a 100644 --- a/arch/univariate/censored_garch.py +++ b/arch/univariate/censored_garch.py @@ -42,11 +42,11 @@ underlying empirical motivation and validation. """ +import math from collections.abc import Sequence from typing import cast import numpy as np -from scipy.stats import norm from arch._typing import ( ArrayLike1D, @@ -56,14 +56,41 @@ ForecastingMethod, RNGType, ) -from arch.univariate.recursions_python import bounds_check_python as bounds_check +from arch.compat.numba import jit +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"] +__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 censored_garch_recursion( + +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, @@ -76,10 +103,15 @@ def censored_garch_recursion( var_bounds: Float64Array2D, ) -> Float64Array1D: """ - Pure-Python reference recursion for :class:`CensoredGARCH`. + Reference recursion for :class:`CensoredGARCH`. - Not currently accelerated (no Cython/Numba backend) -- see the class - docstring's "Notes" section. + 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 ---------- @@ -142,19 +174,24 @@ def censored_garch_recursion( if a <= 0: corr = 1.0 else: - one_minus_phi = norm.sf(a) + 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 + 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""" Censored/latent-variance GARCH(p, q) for return series subject to a @@ -206,13 +243,24 @@ class CensoredGARCH(VolatilityProcess, metaclass=AbstractDocStringInheritor): :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. No accelerated - (Cython/Numba) recursion is provided yet -- see - :func:`censored_garch_recursion`. + 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). + 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. """ def __init__( diff --git a/doc/source/changes/8.0.rst b/doc/source/changes/8.0.rst index 813a898681..b7f5750e35 100644 --- a/doc/source/changes/8.0.rst +++ b/doc/source/changes/8.0.rst @@ -8,8 +8,15 @@ 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). + 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. \ No newline at end of file +- Moved to Meson build system from setuptools. +- Minimum required Python is now 3.10. +- Improved compatability with future NumPy and pandas, including pandas 3. +- Many small documentation fixes. +- Improved typing accuracy across the codebase. + + From 0910e72913bc6911dcc779d1e70acc433e92b647 Mon Sep 17 00:00:00 2001 From: Kevin Sheppard Date: Thu, 8 Oct 2026 06:58:38 +0100 Subject: [PATCH 4/6] Refactor imports and organize dependencies --- arch/univariate/censored_garch.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/arch/univariate/censored_garch.py b/arch/univariate/censored_garch.py index 881163ab5a..c7c0c1d357 100644 --- a/arch/univariate/censored_garch.py +++ b/arch/univariate/censored_garch.py @@ -42,8 +42,10 @@ underlying empirical motivation and validation. """ -import math +from arch.compat.numba import jit + from collections.abc import Sequence +import math from typing import cast import numpy as np @@ -56,11 +58,11 @@ ForecastingMethod, RNGType, ) -from arch.compat.numba import jit 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) From 19eeb578331d846bbb8068f7d0a5429b6f1d5311 Mon Sep 17 00:00:00 2001 From: Muhammad Uzair Date: Thu, 8 Oct 2026 12:43:27 +0500 Subject: [PATCH 5/6] TST: Cover remaining branches in CensoredGARCH Add tests for the length-mismatch error in compute_variance, p=0 / q=0 and higher-order starting values, and the simulate initial-value paths (user-supplied, stationary default, unit-root fallback). Run the deep-tail, zero-threshold and normal-helper tests against the pure-Python sources too, so their lines are counted when numba JIT is active. --- arch/tests/univariate/test_censored_garch.py | 121 ++++++++++++++++--- 1 file changed, 101 insertions(+), 20 deletions(-) diff --git a/arch/tests/univariate/test_censored_garch.py b/arch/tests/univariate/test_censored_garch.py index 7790861a9e..49b0502844 100644 --- a/arch/tests/univariate/test_censored_garch.py +++ b/arch/tests/univariate/test_censored_garch.py @@ -136,7 +136,10 @@ def test_recursion_deep_tail_underflow_branch(): 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.""" - from arch.univariate.censored_garch import censored_garch_recursion + from arch.univariate.censored_garch import ( + censored_garch_recursion, + censored_garch_recursion_python, + ) from scipy.stats import norm a = 40.0 # norm.sf(40) underflows to exactly 0.0 in float64 @@ -149,15 +152,17 @@ def test_recursion_deep_tail_underflow_branch(): resids = np.array([threshold, 0.0]) censored = np.array([True, False]) thresh_arr = np.array([threshold, 0.0]) - sigma2 = np.zeros(2) var_bounds = np.array([[1e-12, 1e12], [1e-12, 1e12]]) - censored_garch_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) + # 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(): @@ -165,7 +170,10 @@ def test_recursion_zero_threshold_branch(): 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.""" - from arch.univariate.censored_garch import censored_garch_recursion + from arch.univariate.censored_garch import ( + censored_garch_recursion, + censored_garch_recursion_python, + ) # omega=0, alpha=1, backcast=1.0 => sigma2[0] = 1.0 exactly (isolates the # step cleanly, same trick as the deep-tail test above). @@ -173,16 +181,18 @@ def test_recursion_zero_threshold_branch(): resids = np.array([0.0, 0.0]) censored = np.array([True, False]) threshold = np.array([0.0, 0.0]) - sigma2 = np.zeros(2) var_bounds = np.array([[1e-12, 1e12], [1e-12, 1e12]]) - censored_garch_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) + # 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(): @@ -211,16 +221,24 @@ def test_norm_helpers_match_scipy(): 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.""" - from arch.univariate.censored_garch import _norm_pdf, _norm_sf + from arch.univariate.censored_garch import ( + _norm_pdf, + _norm_pdf_python, + _norm_sf, + _norm_sf_python, + ) 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: - assert_allclose(_norm_pdf(x), norm.pdf(x), rtol=1e-12, atol=1e-300) + # 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: - assert_allclose(_norm_sf(x), scipy_sf, rtol=1e-8) + for sf in (_norm_sf, _norm_sf_python): + assert_allclose(sf(x), scipy_sf, rtol=1e-8) def test_jitted_recursion_matches_pure_python(): @@ -349,3 +367,66 @@ def test_unsupported_forecast_options(): 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) From c51b360a2d3586824e6d7549005a423b76d3b690 Mon Sep 17 00:00:00 2001 From: Muhammad Uzair Date: Sun, 18 Oct 2026 13:38:04 +0500 Subject: [PATCH 6/6] STY: Fix lint, formatting and numpydoc issues in CensoredGARCH --- arch/tests/univariate/test_censored_garch.py | 48 +++++++++----------- arch/univariate/censored_garch.py | 48 ++++++++++---------- 2 files changed, 46 insertions(+), 50 deletions(-) diff --git a/arch/tests/univariate/test_censored_garch.py b/arch/tests/univariate/test_censored_garch.py index 49b0502844..a9f8df7874 100644 --- a/arch/tests/univariate/test_censored_garch.py +++ b/arch/tests/univariate/test_censored_garch.py @@ -1,11 +1,21 @@ 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 +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") @@ -73,7 +83,6 @@ def test_mismatched_threshold_length(): 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.""" - from arch.univariate.volatility import GARCH resids = setup["resids"] t = setup["t"] @@ -136,11 +145,6 @@ def test_recursion_deep_tail_underflow_branch(): 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.""" - from arch.univariate.censored_garch import ( - censored_garch_recursion, - censored_garch_recursion_python, - ) - from scipy.stats import norm a = 40.0 # norm.sf(40) underflows to exactly 0.0 in float64 assert norm.sf(a) == 0.0 @@ -170,10 +174,6 @@ def test_recursion_zero_threshold_branch(): 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.""" - from arch.univariate.censored_garch import ( - censored_garch_recursion, - censored_garch_recursion_python, - ) # omega=0, alpha=1, backcast=1.0 => sigma2[0] = 1.0 exactly (isolates the # step cleanly, same trick as the deep-tail test above). @@ -200,7 +200,6 @@ def test_recursion_zero_variance_branch(): (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.""" - from arch.univariate.censored_garch import censored_garch_recursion parameters = np.array([0.0, 0.0, 0.0]) # omega=0 -> sigma2[0] forced to 0 resids = np.array([1.5]) @@ -221,12 +220,6 @@ def test_norm_helpers_match_scipy(): 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.""" - from arch.univariate.censored_garch import ( - _norm_pdf, - _norm_pdf_python, - _norm_sf, - _norm_sf_python, - ) 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: @@ -246,10 +239,6 @@ def test_jitted_recursion_matches_pure_python(): 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.""" - from arch.univariate.censored_garch import ( - censored_garch_recursion, - censored_garch_recursion_python, - ) rng = RandomState(7) nobs = 2000 @@ -264,7 +253,16 @@ def test_jitted_recursion_matches_pure_python(): 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 + 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 @@ -277,7 +275,6 @@ def test_correction_matches_truncated_normal_identity(): 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.""" - from arch.univariate.censored_garch import censored_garch_recursion omega = 1.0 threshold = 1.2 @@ -315,7 +312,6 @@ def test_correction_matches_truncated_normal_identity(): 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.""" - import pandas as pd rng = RandomState(42) nobs = 500 @@ -382,7 +378,7 @@ def test_compute_variance_length_mismatch(): ) -@pytest.mark.parametrize("p, q", [(1, 0), (0, 1), (2, 2)]) +@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.""" diff --git a/arch/univariate/censored_garch.py b/arch/univariate/censored_garch.py index c7c0c1d357..df4494095e 100644 --- a/arch/univariate/censored_garch.py +++ b/arch/univariate/censored_garch.py @@ -62,8 +62,11 @@ 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"] +__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) @@ -118,7 +121,7 @@ def censored_garch_recursion_python( Parameters ---------- parameters : ndarray - [omega, alpha_1..alpha_p, beta_1..beta_q] + 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 @@ -196,8 +199,10 @@ def censored_garch_recursion_python( class CensoredGARCH(VolatilityProcess, metaclass=AbstractDocStringInheritor): r""" - Censored/latent-variance GARCH(p, q) for return series subject to a - known, time-varying censoring rule (e.g. an exchange price-limit band). + 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 ---------- @@ -215,12 +220,6 @@ class CensoredGARCH(VolatilityProcess, metaclass=AbstractDocStringInheritor): q : int, optional Order of the lagged conditional variance. Default is 1. - 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) - Notes ----- In this class of processes, the variance dynamics are @@ -263,6 +262,12 @@ class CensoredGARCH(VolatilityProcess, metaclass=AbstractDocStringInheritor): 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__( @@ -402,14 +407,11 @@ def simulate( burn: int = 500, initial_value: float | Float64Array | None = None, ) -> tuple[Float64Array, Float64Array]: - """ - Simulate from the model, assuming NO censoring in the simulated - path (future censoring status is not known ex ante -- this - simulates the *unconditional* process the model implies, which is - only exactly correct when the true return process is never - censored; it is an approximation whenever a future censoring rule - would actually bind). See class Notes. - """ + # 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) @@ -451,8 +453,7 @@ def _check_forecasting_method( ) if method != "analytic": raise NotImplementedError( - "Only method='analytic' is currently supported for " - "CensoredGARCH." + "Only method='analytic' is currently supported for CensoredGARCH." ) def _analytic_forecast( @@ -464,7 +465,7 @@ def _analytic_forecast( start: int, horizon: int, ) -> VarianceForecast: - sigma2, forecasts = self._one_step_forecast( + _, forecasts = self._one_step_forecast( parameters, to_array_1d(resids), backcast, var_bounds, horizon, start ) return VarianceForecast(forecasts) @@ -481,6 +482,5 @@ def _simulation_forecast( rng: RNGType, ) -> VarianceForecast: raise NotImplementedError( - "Simulation-based forecasts are not yet implemented for " - "CensoredGARCH." + "Simulation-based forecasts are not yet implemented for CensoredGARCH." )