BUG: Use the expected variance path in RiskMetrics2006 forecasts

arch/tests/univariate/test_variance_forecasting.py

@@ -1808,8 +1808,20 @@ def test_rm2006_forecast(self):
         assert forecasts.forecast_paths is None
         assert forecasts.shocks is None
         one_step = _sigma2[1:]
+        assert_allclose(forecasts.forecasts[:, 0], one_step, rtol=1e-4)
+        # The components are EWMAs of squared residuals, and the expected
+        # squared residual is the combined variance
+        w = vol._ewma_combination_weights()
+        mus = vol._ewma_smoothing_parameters()
+        resids2 = np.asarray(resids) ** 2
+        components = np.empty((resids2.shape[0] + 1, mus.shape[0]))
+        components[0] = backcast
+        for t in range(resids2.shape[0]):
+            components[t + 1] = mus * components[t] + (1 - mus) * resids2[t]
+        expected = components[1:]
         for i in range(10):
-            assert_allclose(forecasts.forecasts[:, i], one_step, rtol=1e-4)
+            assert_allclose(forecasts.forecasts[:, i], expected @ w, rtol=1e-4)
+            expected = mus * expected + (1 - mus) * (expected @ w)[:, None]
 
         alt_forecasts = vol.forecast(
             params, resids, backcast, var_bounds, horizon=10, start=500

arch/univariate/volatility.py

@@ -2472,8 +2472,34 @@ def _analytic_forecast(
             horizon,
             start_index=start,
         )
+        if horizon == 1:
+            return VarianceForecast(forecasts)
+
+        kmax = self.kmax
+        w = self._ewma_combination_weights()
+        mus = self._ewma_smoothing_parameters()
+        backcast = cast("Float64Array1D", to_array_1d(np.asarray(backcast)))
+        t = resids.shape[0]
+        components = np.empty((kmax, t + 1))
+        _resids = np.ascontiguousarray(resids)
+        for k in range(kmax):
+            ewma_recursion(
+                mus[k],
+                to_array_1d(_resids),
+                to_array_1d(components[k, :]),
+                t + 1,
+                backcast[k],
+            )
+        # Each EWMA component is updated with the squared residual, whose
+        # expected value is the combined variance, so the expected values of
+        # the components follow a linear recursion. The forecasts are then
+        # reweighted combinations of the final component values, with weights
+        # w(h + 1) = mus * w(h) + w * sum((1 - mus) * w(h)) and w(1) = w.
+        components = components.T[start + 1 :]
+        weights = w
         for i in range(1, horizon):
-            forecasts[:, i] = forecasts[:, 0]
+            weights = mus * weights + w * ((1 - mus) @ weights)
+            forecasts[:, i] = components @ weights
         return VarianceForecast(forecasts)
 
     def _simulation_forecast(