ENH: Add FIAPARCH volatility model

arch/tests/univariate/test_forecast.py

@@ -14,6 +14,7 @@
     APARCH,
     ARX,
     EGARCH,
+    FIAPARCH,
     FIGARCH,
     GARCH,
     HARCH,
@@ -42,6 +43,7 @@
     ConstantVariance(),
     GARCH(),
     FIGARCH(),
+    FIAPARCH(),
     EWMAVariance(lam=0.94),
     MIDASHyperbolic(),
     HARCH(lags=[1, 5, 22]),

arch/tests/univariate/test_mean.py

@@ -46,6 +46,7 @@
     APARCH,
     ARCH,
     EGARCH,
+    FIAPARCH,
     FIGARCH,
     GARCH,
     HARCH,
@@ -685,6 +686,9 @@ def test_arch_model(self):
         am = arch_model(self.y, vol="figarch")
         assert isinstance(am.volatility, FIGARCH)
 
+        am = arch_model(self.y, vol="fiaparch")
+        assert isinstance(am.volatility, FIAPARCH)
+
         am = arch_model(self.y, vol="aparch")
         assert isinstance(am.volatility, APARCH)
 
@@ -1354,6 +1358,23 @@ def test_invalid_vol_dist():
         ConstantMean(SP500, distribution="Skew-t")
 
 
+def test_fiaparch_non_int_p():
+    with pytest.raises(TypeError, match=r"p must be a scalar int"):
+        arch_model(SP500, vol="fiaparch", p=[1, 2])
+
+
+def test_fiaparch_arch_model_kwargs():
+    am = arch_model(SP500, vol="fiaparch", o=0)
+    assert isinstance(am.volatility, FIAPARCH)
+    assert am.volatility.o == 0
+    assert am.volatility.name == "FI Power ARCH"
+
+    am2 = arch_model(SP500, vol="FIAPARCH", p=0, q=0)
+    assert isinstance(am2.volatility, FIAPARCH)
+    assert am2.volatility.p == 0
+    assert am2.volatility.q == 0
+
+
 def test_param_cov():
     mod = ConstantMean(SP500)
     res = mod.fit(disp="off")

arch/tests/univariate/test_variance_forecasting.py

@@ -16,6 +16,7 @@
 from arch.univariate.volatility import (
     APARCH,
     EGARCH,
+    FIAPARCH,
     FIGARCH,
     GARCH,
     HARCH,
@@ -1930,6 +1931,141 @@ def test_aparch_simulation_smoke(self, o, delta):
                 method="analytic",
             )
 
+    def test_fiaparch_one_step(self):
+        trunc = 50
+        vol = FIAPARCH(truncation=trunc)
+        resids = self.resid
+        backcast = vol.backcast(resids)
+        var_bounds = vol.variance_bounds(resids)
+        params = np.array([0.1, 0.2, 0.4, 0.2, -0.3, 1.5])
+        sigma2 = np.empty_like(resids)
+        vol.compute_variance(params, resids, sigma2, backcast, var_bounds)
+        forecast = vol.forecast(
+            params, resids, backcast, var_bounds, horizon=1, start=0
+        )
+        assert_allclose(sigma2[1:], forecast.forecasts[:-1, 0])
+
+        delta = 1.5
+        vol_fixed = FIAPARCH(truncation=trunc, delta=delta)
+        params_fixed = np.array([0.1, 0.2, 0.4, 0.2, -0.3])
+        sigma2_f = np.empty_like(resids)
+        vol_fixed.compute_variance(
+            params_fixed, resids, sigma2_f, backcast, var_bounds
+        )
+        forecast_f = vol_fixed.forecast(
+            params_fixed, resids, backcast, var_bounds, horizon=1, start=0
+        )
+        assert_allclose(sigma2_f[1:], forecast_f.forecasts[:-1, 0])
+
+        with pytest.raises(
+            ValueError, match=r"Analytic forecasts not available for horizon"
+        ):
+            vol.forecast(
+                params, resids, backcast, var_bounds, horizon=2, method="analytic"
+            )
+
+    @pytest.mark.parametrize("o", [0, 1])
+    @pytest.mark.parametrize("delta", [None, 1.5])
+    def test_fiaparch_simulation_smoke(self, o, delta):
+        dist = Normal(seed=self.rng)
+        rng = dist.simulate([])
+        trunc = 50
+        vol = FIAPARCH(o=o, delta=delta, truncation=trunc)
+        resids = self.resid
+        backcast = vol.backcast(resids)
+        var_bounds = vol.variance_bounds(resids)
+        params = [0.1, 0.2, 0.4, 0.2]
+        if o == 0:
+            params = [0.1, 0.2, 0.4, 0.2]
+        else:
+            params = [0.1, 0.2, 0.4, 0.2, -0.3]
+        if delta is None:
+            params = np.array(params + [1.5])
+        else:
+            params = np.array(params)
+        sigma2 = np.empty_like(resids)
+        vol.compute_variance(params, resids, sigma2, backcast, var_bounds)
+        forecast = vol.forecast(
+            params,
+            resids,
+            backcast,
+            var_bounds,
+            horizon=10,
+            start=0,
+            method="simulation",
+            rng=rng,
+            simulations=100,
+        )
+        assert_allclose(sigma2[1:], forecast.forecasts[:-1, 0])
+        assert forecast.forecast_paths is not None
+        assert forecast.shocks is not None
+        with pytest.raises(ValueError, match=r"Analytic forecasts not"):
+            vol.forecast(
+                params,
+                resids,
+                backcast,
+                var_bounds,
+                horizon=10,
+                start=0,
+                method="analytic",
+            )
+
+    @pytest.mark.parametrize("p,q", [(0, 1), (1, 0), (0, 0)])
+    def test_fiaparch_one_step_reduced(self, p, q):
+        trunc = 50
+        vol = FIAPARCH(p=p, q=q, truncation=trunc)
+        resids = self.resid
+        backcast = vol.backcast(resids)
+        var_bounds = vol.variance_bounds(resids)
+        params = [0.1]
+        if p:
+            params.append(0.2)
+        params.append(0.4)
+        if q:
+            params.append(0.2)
+        params.extend([-0.3, 1.5])
+        params = np.array(params)
+        sigma2 = np.empty_like(resids)
+        vol.compute_variance(params, resids, sigma2, backcast, var_bounds)
+        forecast = vol.forecast(
+            params, resids, backcast, var_bounds, horizon=1, start=0
+        )
+        assert_allclose(sigma2[1:], forecast.forecasts[:-1, 0])
+
+    @pytest.mark.parametrize("p,q", [(0, 1), (1, 0), (0, 0)])
+    def test_fiaparch_simulation_reduced(self, p, q):
+        dist = Normal(seed=self.rng)
+        rng = dist.simulate([])
+        trunc = 50
+        vol = FIAPARCH(p=p, q=q, truncation=trunc)
+        resids = self.resid
+        backcast = vol.backcast(resids)
+        var_bounds = vol.variance_bounds(resids)
+        params = [0.1]
+        if p:
+            params.append(0.2)
+        params.append(0.4)
+        if q:
+            params.append(0.2)
+        params.extend([-0.3, 1.5])
+        params = np.array(params)
+        sigma2 = np.empty_like(resids)
+        vol.compute_variance(params, resids, sigma2, backcast, var_bounds)
+        forecast = vol.forecast(
+            params,
+            resids,
+            backcast,
+            var_bounds,
+            horizon=10,
+            start=0,
+            method="simulation",
+            rng=rng,
+            simulations=100,
+        )
+        assert_allclose(sigma2[1:], forecast.forecasts[:-1, 0])
+        assert forecast.forecast_paths is not None
+        assert forecast.shocks is not None
+
     def test_midas_analytical(self):
         vol = MIDASHyperbolic()
         resids = self.resid