BUG: Fix bugs in fixed result

Systematically find an fix bugs that affect fixed results

closes #749

Author: bashtage

arch/tests/univariate/test_mean.py

@@ -1,6 +1,7 @@
 from arch.compat.pandas import MONTH_END
 
 from io import StringIO
+import itertools
 from itertools import product
 from string import ascii_lowercase
 import struct
@@ -26,7 +27,12 @@
 
 from arch.data import sp500
 from arch.typing import Literal
-from arch.univariate.base import ARCHModelForecast, ARCHModelResult, _align_forecast
+from arch.univariate.base import (
+    ARCHModel,
+    ARCHModelForecast,
+    ARCHModelResult,
+    _align_forecast,
+)
 from arch.univariate.distribution import (
     GeneralizedError,
     Normal,
@@ -46,6 +52,7 @@
     FixedVariance,
     MIDASHyperbolic,
     RiskMetrics2006,
+    VolatilityProcess,
 )
 from arch.utility.exceptions import ConvergenceWarning, DataScaleWarning
 
@@ -73,6 +80,14 @@
 DISPLAY: Literal["off", "final"] = "off"
 UPDATE_FREQ = 0 if DISPLAY == "off" else 3
 SP500 = 100 * sp500.load()["Adj Close"].pct_change().dropna()
+rs = np.random.RandomState(20241029)
+X = SP500 * 0.01 + SP500.std() * rs.standard_normal(SP500.shape)
+
+
+def close_plots():
+    import matplotlib.pyplot as plt
+
+    plt.close("all")
 
 
 @pytest.fixture(scope="module", params=[True, False])
@@ -83,6 +98,73 @@ def simulated_data(request):
     return np.asarray(sim_data.data) if request.param else sim_data.data
 
 
+simple_mean_models = [
+    ARX(SP500, lags=1),
+    HARX(SP500, lags=[1, 5]),
+    ConstantMean(SP500),
+    ZeroMean(SP500),
+]
+
+mean_models = [
+    ARX(SP500, x=X, lags=1),
+    HARX(SP500, x=X, lags=[1, 5]),
+    LS(SP500, X),
+] + simple_mean_models
+
+analytic_volatility_processes = [
+    ARCH(3),
+    FIGARCH(1, 1),
+    GARCH(1, 1, 1),
+    HARCH([1, 5, 22]),
+    ConstantVariance(),
+    EWMAVariance(0.94),
+    FixedVariance(np.full_like(SP500, SP500.var())),
+    MIDASHyperbolic(),
+    RiskMetrics2006(),
+]
+
+other_volatility_processes = [
+    APARCH(1, 1, 1, 1.5),
+    EGARCH(1, 1, 1),
+]
+
+volatility_processes = analytic_volatility_processes + other_volatility_processes
+
+
+@pytest.fixture(
+    scope="module",
+    params=list(itertools.product(simple_mean_models, analytic_volatility_processes)),
+    ids=[
+        f"{a.__class__.__name__}-{b}"
+        for a, b in itertools.product(simple_mean_models, analytic_volatility_processes)
+    ],
+)
+def forecastable_model(request):
+    mod: ARCHModel
+    vol: VolatilityProcess
+    mod, vol = request.param
+    mod.volatility = vol
+    res = mod.fit()
+    return res, mod.fix(res.params)
+
+
+@pytest.fixture(
+    scope="module",
+    params=list(itertools.product(mean_models, volatility_processes)),
+    ids=[
+        f"{a.__class__.__name__}-{b}"
+        for a, b in itertools.product(mean_models, volatility_processes)
+    ],
+)
+def fit_fixed_models(request):
+    mod: ARCHModel
+    vol: VolatilityProcess
+    mod, vol = request.param
+    mod.volatility = vol
+    res = mod.fit()
+    return res, mod.fix(res.params)
+
+
 class TestMeanModel:
     @classmethod
     def setup_class(cls):
@@ -480,9 +562,7 @@ def test_ar_plot(self):
         with pytest.raises(ValueError):
             res.plot(annualize="unknown")
 
-        import matplotlib.pyplot as plt
-
-        plt.close("all")
+        close_plots()
 
         res.plot(scale=360)
         res.hedgehog_plot(start=500)
@@ -491,7 +571,7 @@ def test_ar_plot(self):
         res.hedgehog_plot(start=500, method="simulation", simulations=100)
         res.hedgehog_plot(plot_type="volatility", method="bootstrap")
 
-        plt.close("all")
+        close_plots()
 
     def test_arch_arx(self):
         self.rng.seed(12345)
@@ -1370,3 +1450,61 @@ def test_non_contiguous_input(use_numpy):
     mod = arch_model(y, mean="Zero")
     res = mod.fit()
     assert res.params.shape[0] == 3
+
+
+def test_fixed_equivalence(fit_fixed_models):
+    res, res_fixed = fit_fixed_models
+
+    assert_allclose(res.aic, res_fixed.aic)
+    assert_allclose(res.bic, res_fixed.bic)
+    assert_allclose(res.loglikelihood, res_fixed.loglikelihood)
+    assert res.nobs == res_fixed.nobs
+    assert res.num_params == res_fixed.num_params
+    assert_allclose(res.params, res_fixed.params)
+    assert_allclose(res.conditional_volatility, res_fixed.conditional_volatility)
+    assert_allclose(res.std_resid, res_fixed.std_resid)
+    assert_allclose(res.resid, res_fixed.resid)
+    assert_allclose(res.arch_lm_test(5).stat, res_fixed.arch_lm_test(5).stat)
+    assert res.model.__class__ is res_fixed.model.__class__
+    assert res.model.volatility.__class__ is res_fixed.model.volatility.__class__
+    assert isinstance(res.summary(), type(res_fixed.summary()))
+    if res.num_params > 0:
+        assert "std err" in str(res.summary())
+        assert "std err" not in str(res_fixed.summary())
+
+
+@pytest.mark.skipif(not HAS_MATPLOTLIB, reason="matplotlib not installed")
+def test_fixed_equivalence_plots(fit_fixed_models):
+    res, res_fixed = fit_fixed_models
+
+    fig = res.plot()
+    fixed_fig = res_fixed.plot()
+    assert isinstance(fig, type(fixed_fig))
+
+    close_plots()
+
+
+@pytest.mark.parametrize("simulations", [1, 100])
+def test_fixed_equivalence_forecastable(forecastable_model, simulations):
+    res, res_fixed = forecastable_model
+    f1 = res.forecast(horizon=5)
+    f2 = res_fixed.forecast(horizon=5)
+    assert isinstance(f1, type(f2))
+    assert_allclose(f1.mean, f2.mean)
+    assert_allclose(f1.variance, f2.variance)
+
+    f1 = res.forecast(horizon=5, method="simulation", simulations=simulations)
+    f2 = res_fixed.forecast(horizon=5, method="simulation", simulations=simulations)
+    assert isinstance(f1, type(f2))
+    f1 = res.forecast(horizon=5, method="bootstrap", simulations=simulations)
+    f2 = res_fixed.forecast(horizon=5, method="bootstrap", simulations=simulations)
+    assert isinstance(f1, type(f2))
+
+
+@pytest.mark.skipif(not HAS_MATPLOTLIB, reason="matplotlib not installed")
+def test_fixed_equivalence_forecastable_plots(forecastable_model):
+    res, res_fixed = forecastable_model
+    fig1 = res.hedgehog_plot(start=SP500.shape[0] - 25)
+    fig2 = res_fixed.hedgehog_plot(start=SP500.shape[0] - 25)
+    assert isinstance(fig1, type(fig2))
+    close_plots()

arch/univariate/base.py

@@ -386,9 +386,11 @@ def _fit_parameterless_model(
         vol_final[first_obs:last_obs] = vol
 
         names = self._all_parameter_names()
-        loglikelihood = self._static_gaussian_loglikelihood(y)
         r2 = self._r2(params)
         fit_start, fit_stop = self._fit_indices
+        loglikelihood = -1.0 * self._loglikelihood(
+            params, vol**2 * np.ones(fit_stop - fit_start), backcast, var_bounds
+        )
 
         assert isinstance(r2, float)
         return ARCHModelResult(
@@ -523,8 +525,7 @@ def fix(
         names = self._all_parameter_names()
         # Reshape resids and vol
         first_obs, last_obs = self._fit_indices
-        resids_final = np.empty_like(self._y, dtype=np.double)
-        resids_final.fill(np.nan)
+        resids_final = np.full_like(self._y, np.nan, dtype=np.double)
         resids_final[first_obs:last_obs] = resids
         vol_final = np.empty_like(self._y, dtype=np.double)
         vol_final.fill(np.nan)
@@ -533,8 +534,8 @@ def fix(
         model_copy = deepcopy(self)
         return ARCHModelFixedResult(
             params,
-            resids,
-            vol,
+            resids_final,
+            vol_final,
             self._y_series,
             names,
             loglikelihood,

arch/univariate/mean.py

@@ -1018,7 +1018,11 @@ def forecast(
             for i in range(horizon):
                 _impulses = impulse[i::-1][:, None]
                 lrvp = variance_paths[:, :, : (i + 1)].dot(_impulses**2)
-                long_run_variance_paths[:, :, i] = np.squeeze(lrvp)
+                lrvp = np.squeeze(lrvp)
+                if lrvp.ndim < 2:
+                    lrvp = np.atleast_1d(lrvp)
+                    lrvp = lrvp[None, :]
+                long_run_variance_paths[:, :, i] = lrvp
             t, m = self._y.shape[0], self._max_lags
             mean_paths = np.empty(shocks.shape[:2] + (m + horizon,))
             dynp_rev = dynp[::-1]

arch/univariate/volatility.py

@@ -930,8 +930,7 @@ def _analytic_forecast(
         forecasts = np.full((t - start, horizon), np.nan)
 
         forecasts[:, :] = parameters[0]
-        forecast_paths = None
-        return VarianceForecast(forecasts, forecast_paths)
+        return VarianceForecast(forecasts)
 
     def _simulation_forecast(
         self,
@@ -2804,7 +2803,8 @@ class FixedVariance(VolatilityProcess, metaclass=AbstractDocStringInheritor):
     ----------
     variance : {array, Series}
         Array containing the variances to use.  Should have the same shape as the
-        data used in the model.
+        data used in the model. This is not checked since the model is not
+        available when the FixedVariance process is created.
     unit_scale : bool, optional
         Flag whether to enforce a unit scale.  If False, a scale parameter will be
         estimated so that the model variance will be proportional to ``variance``.
@@ -2823,7 +2823,7 @@ def __init__(self, variance: Float64Array, unit_scale: bool = False) -> None:
         self._name = "Fixed Variance"
         self._name += " (Unit Scale)" if unit_scale else ""
         self._variance_series = ensure1d(variance, "variance", True)
-        self._variance = np.asarray(variance)
+        self._variance = np.atleast_1d(variance)
 
     def compute_variance(
         self,
@@ -2841,6 +2841,10 @@ def compute_variance(
         return sigma2
 
     def starting_values(self, resids: Float64Array) -> Float64Array:
+        if self._variance.ndim != 1 or self._variance.shape[0] < self._stop:
+            raise ValueError(
+                f"variance must be a 1-d array with at least {self._stop} elements"
+            )
         if not self._unit_scale:
             _resids = resids / np.sqrt(self._variance[self._start : self._stop])
             return np.array([_resids.var()])
@@ -2893,7 +2897,7 @@ def _analytic_forecast(
         horizon: int,
     ) -> VarianceForecast:
         t = resids.shape[0]
-        forecasts = np.full((t, horizon), np.nan)
+        forecasts = np.full((t - start, horizon), np.nan)
 
         return VarianceForecast(forecasts)
 
@@ -2909,10 +2913,9 @@ def _simulation_forecast(
         rng: RNGType,
     ) -> VarianceForecast:
         t = resids.shape[0]
-        forecasts = np.full((t, horizon), np.nan)
-        forecast_paths = np.empty((t, simulations, horizon))
-        forecast_paths.fill(np.nan)
-        shocks = np.full((t, simulations, horizon), np.nan)
+        forecasts = np.full((t - start, horizon), np.nan)
+        forecast_paths = np.full((t - start, simulations, horizon), np.nan)
+        shocks = np.full((t - start, simulations, horizon), np.nan)
 
         return VarianceForecast(forecasts, forecast_paths, shocks)