add polymorphic compute_filtered_probs to VolatilityProcess and MSGARCH. store filtered_probs in ARCHModelResult

Author: wday0507

arch/tests/univariate/test_volatility.py

@@ -10,6 +10,7 @@
     assert_array_equal,
     assert_equal,
 )
+from sklearn.mixture import GaussianMixture
 import pytest
 from scipy.special import gamma, gammaln
 from arch import arch_model

arch/univariate/base.py

@@ -863,10 +863,7 @@ def fit(
         resids_final = np.full(self._y.shape, np.nan, dtype=float)
         resids_final[first_obs:last_obs] = resids
 
-        if isinstance(self.volatility, MSGARCH):
-            filtered_probs = self.volatility.filtered_probs  # n_regimes x n_obs
-        else:
-            filtered_probs = None
+        filtered_probs = self.volatility.compute_filtered_probs(params, resids, sigma2, backcast, var_bounds)
 
 
         if isinstance(self.volatility, MSGARCH):   
@@ -896,6 +893,7 @@ def fit(
             fit_start,
             fit_stop,
             model_copy,
+            filtered_probs,
         )
 
     @abstractmethod
@@ -1829,6 +1827,7 @@ def __init__(
         fit_start: int,
         fit_stop: int,
         model: ARCHModel,
+        filtered_probs: Float64Array | None = None,
     ) -> None:
         super().__init__(
             params, resid, volatility, dep_var, names, loglikelihood, is_pandas, model
@@ -1839,6 +1838,7 @@ def __init__(
         self._r2 = r2
         self.cov_type: str = cov_type
         self._optim_output = optim_output
+        self._filtered_probs = filtered_probs
 
     @cached_property
     def scale(self) -> float:
@@ -2089,6 +2089,10 @@ def optimization_result(self) -> OptimizeResult:
             Result from numerical optimization of the log-likelihood.
         """
         return self._optim_output
+    
+    @property
+    def filtered_probs(self):
+        return self._filtered_probs 
 
 
 def _align_forecast(

arch/univariate/volatility.py

@@ -861,6 +861,12 @@ def parameter_names(self) -> list[str]:
         """
 
 
+    def compute_filtered_probs(self, params, resids, sigma2, backcast, var_bounds):
+        """
+        Default behavior for non MS models. Returns 1's for each time period (we are always in the single regime).
+        """
+        return np.ones((1, len(resids)))
+
 class ConstantVariance(VolatilityProcess, metaclass=AbstractDocStringInheritor):
     r"""
     Constant volatility process
@@ -4011,10 +4017,6 @@ def msgarch_loglikelihood(self, parameters, resids, sigma2, backcast, var_bounds
         #mean params, volatility params per regime, distribution params
         mp, vp, dp = self._parse_parameters(parameters)
 
-        # GARCH params per regime
-        garch_block = vp[:self.k*self.num_params_single]
-        garch_params = garch_block.reshape(self.k, self.num_params_single)  
-
         # extract the free parameters for the transition matrix
         start = self.k * self.num_params_single
         end = start + self.k  # only one free param per row
@@ -4110,5 +4112,46 @@ def _parse_parameters(self, x: np.ndarray):
         dp = x[km+kv:]
 
         return mp, vp, dp
+    
+
+    def compute_filtered_probs(self, parameters, resids, sigma2, backcast, var_bounds):
+        # parse parameters
+        mp, vp, dp = self._parse_parameters(parameters)
+
+        # transition matrix
+        start = self.k * self.num_params_single
+        end = start + self.k
+        free_P = vp[start:end]
+        transition_matrix = np.zeros((self.k, self.k))
+        transition_matrix[0, 0] = free_P[0]
+        transition_matrix[0, 1] = 1.0 - free_P[0]
+        transition_matrix[1, 0] = free_P[1]
+        transition_matrix[1, 1] = 1.0 - free_P[1]
+
+        #  Initial probabilties
+        eigvals, eigvecs = np.linalg.eig(transition_matrix.T)
+        pi = np.real(eigvecs[:, np.isclose(eigvals, 1)].flatten())
+        pi /= pi.sum()
+
+        
+        p = np.zeros((self.k, len(resids)))
+        p[:, 0] = pi
+        sigma2 = self.compute_variance(vp, resids, np.zeros((len(resids), self.k)), backcast, var_bounds)
+
+        # Hamilton filter
+        var_floor = 1e-4
+        sigma2 = np.maximum(sigma2, var_floor)
+        log_p = np.log(np.maximum(p, 1e-12))
+
+        for t in range(1, len(resids)):
+            sigma2_t = sigma2[t, :]
+            log_likes = -0.5 * (np.log(2 * np.pi * sigma2_t) + (resids[t] ** 2) / sigma2_t)
+            log_pred_probs = logsumexp(np.log(np.maximum(transition_matrix.T, 1e-12)) + log_p[:, t-1][:, None], axis=0)
+            log_mixture = logsumexp(log_pred_probs + log_likes)
+            log_p[:, t] = log_pred_probs + log_likes - log_mixture
+
+        # return filtered probs
+        return np.exp(log_p)
+