Add Markov Switching GARCH volatility model

arch/tests/univariate/test_volatility.py

@@ -10,9 +10,10 @@
     assert_array_equal,
     assert_equal,
 )
+from sklearn.mixture import GaussianMixture
 import pytest
 from scipy.special import gamma, gammaln
-
+from arch import arch_model
 from arch.univariate import ZeroMean
 from arch.univariate.distribution import Normal, SkewStudent, StudentsT
 import arch.univariate.recursions_python as recpy
@@ -28,6 +29,7 @@
     FixedVariance,
     MIDASHyperbolic,
     RiskMetrics2006,
+    MSGARCH,
 )
 from arch.utility.exceptions import InitialValueWarning
 
@@ -1804,3 +1806,28 @@ def test_figarch_weights():
     lam_cy = rec.figarch_weights(params, 1, 1, 1000)
     assert_allclose(lam_py, lam_nb)
     assert_allclose(lam_py, lam_cy)
+
+
+
+def test_msgarch():
+    data = np.random.randn(100) # arbitary data series
+    model = arch_model(data, vol="MSGARCH",p=1,q=1,mean="zero")  # 2 regimes
+    result = model.fit(disp="off")
+    forecast = result.forecast(horizon=1)
+    cond_var = result.conditional_volatility
+    probs = result.model.volatility.filtered_probs
+    ll = result.loglikelihood
+
+    assert hasattr(result, "params")
+    assert result is not None
+    assert result.params.shape[0] == 8 #  2*3 GARCH params + 2 transition matrix probabilities (this is valid whilst k=2)
+    assert np.all(np.isfinite(result.params))
+    assert result.model.volatility.k == 2 
+    assert np.all(forecast.variance.values > 0)
+    assert forecast.variance.shape[1] == 1
+    assert np.allclose(probs[:,1:].sum(axis=0), 1.0, atol=1e-6) # excludes t=0 prob's as these have not been filtered
+    assert np.all((probs >= 0) & (probs <= 1))
+    assert cond_var.shape[0] == len(data)
+    assert np.all(cond_var > 0)
+    assert np.isfinite(ll)
+

arch/univariate/__init__.py

@@ -29,6 +29,7 @@
     FixedVariance,
     MIDASHyperbolic,
     RiskMetrics2006,
+    MSGARCH,
 )
 
 recursions: types.ModuleType
@@ -55,6 +56,7 @@
     "HARX",
     "LS",
     "MIDASHyperbolic",
+    "MSGARCH",
     "Normal",
     "RiskMetrics2006",
     "SkewStudent",

arch/univariate/base.py

@@ -34,7 +34,7 @@
     Literal,
 )
 from arch.univariate.distribution import Distribution, Normal
-from arch.univariate.volatility import ConstantVariance, VolatilityProcess
+from arch.univariate.volatility import ConstantVariance, VolatilityProcess, MSGARCH
 from arch.utility.array import ensure1d, to_array_1d
 from arch.utility.exceptions import (
     ConvergenceWarning,
@@ -714,15 +714,21 @@ def fit(
         if total_params == 0:
             return self._fit_parameterless_model(cov_type=cov_type, backcast=backcast)
 
-        sigma2 = np.zeros(resids.shape[0], dtype=float)
-        self._backcast = backcast
-        sv_volatility = v.starting_values(resids)
+        n_regimes = v.k if isinstance(v, MSGARCH) else 1
+        sigma2 = np.zeros((resids.shape[0], n_regimes)) if n_regimes > 1 else np.zeros(resids.shape[0])
+        self._backcast = backcast 
+        sv_volatility = v.starting_values(resids) # initial guess for GARCH recursion
         self._var_bounds = var_bounds = v.variance_bounds(resids)
-        v.compute_variance(sv_volatility, resids, sigma2, backcast, var_bounds)
-        std_resids = resids / np.sqrt(sigma2)
+        v.compute_variance(sv_volatility, resids, sigma2, backcast, var_bounds) # fills sigma 2 recursively using chosen vol model 
+        if n_regimes == 1:
+            std_resids = resids / np.sqrt(sigma2) 
+        else:
+            pi = self.volatility.pi # shape (k,)
+            pi_weighted_sigma2 = sigma2 @ pi # (t,k) @ (k,) = (t,)
+            std_resids = resids / np.sqrt(pi_weighted_sigma2) ## Using stationary distribution to weight regime variances. This is only an approximation (true weights should come from filtered probabilties), but we don't have these available at this stage
 
         # 2. Construct constraint matrices from all models and distribution
-        constraints = (
+        constraints = (                                                               
             self.constraints(),
             self.volatility.constraints(),
             self.distribution.constraints(),
@@ -790,12 +796,16 @@ def fit(
             _callback_info["display"] = update_freq
         disp_flag = True if disp == "final" else False
 
-        func = self._loglikelihood
-        args = (sigma2, backcast, var_bounds)
-        ineq_constraints = constraint(a, b)
+        if isinstance(self.volatility, MSGARCH):  
+            func = self.volatility.msgarch_loglikelihood # MS GARCH
+            args = (resids, sigma2, backcast, var_bounds)
 
-        from scipy.optimize import minimize
+        else:
+            func = self._loglikelihood  # standard GARCH
+            args = (sigma2, backcast, var_bounds)
 
+        ineq_constraints = constraint(a, b)
+        from scipy.optimize import minimize
         options = {} if options is None else options
         options.setdefault("disp", disp_flag)
         with warnings.catch_warnings():
@@ -835,7 +845,7 @@ def fit(
         mp, vp, dp = self._parse_parameters(params)
 
         resids = self.resids(mp)
-        vol = np.zeros(resids.shape[0], dtype=float)
+        vol = self.volatility._initialise_vol(resids, n_regimes)
         self.volatility.compute_variance(vp, resids, vol, backcast, var_bounds)
         vol = cast(Float64Array1D, np.sqrt(vol))
 
@@ -849,8 +859,18 @@ def fit(
         first_obs, last_obs = self._fit_indices
         resids_final = np.full(self._y.shape, np.nan, dtype=float)
         resids_final[first_obs:last_obs] = resids
-        vol_final = np.full(self._y.shape, np.nan, dtype=float)
-        vol_final[first_obs:last_obs] = vol
+
+        filtered_probs = self.volatility.compute_filtered_probs(params, resids, sigma2, backcast, var_bounds)
+
+
+        if isinstance(self.volatility, MSGARCH):   
+            vol_final = np.full(self._y.shape, np.nan, dtype=float)
+            weighted_sigma2 = (vol**2 * filtered_probs.T).sum(axis=1) 
+            vol_final[first_obs:last_obs] = np.sqrt(weighted_sigma2)
+
+        else:              
+            vol_final = np.full(self._y.shape, np.nan, dtype=float)
+            vol_final[first_obs:last_obs] = vol
 
         fit_start, fit_stop = self._fit_indices
         model_copy = deepcopy(self)
@@ -870,6 +890,7 @@ def fit(
             fit_start,
             fit_stop,
             model_copy,
+            filtered_probs,
         )
 
     @abstractmethod
@@ -1803,6 +1824,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
@@ -1813,6 +1835,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:
@@ -2063,6 +2086,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/mean.py

@@ -67,6 +67,7 @@
     HARCH,
     ConstantVariance,
     VolatilityProcess,
+    MSGARCH
 )
 from arch.utility.array import (
     AbstractDocStringInheritor,
@@ -1891,7 +1892,7 @@ def arch_model(
     ] = "Constant",
     lags: int | list[int] | Int32Array | Int64Array | None = 0,
     vol: Literal[
-        "GARCH", "ARCH", "EGARCH", "FIGARCH", "APARCH", "HARCH", "FIGARCH"
+        "GARCH", "ARCH", "EGARCH", "FIGARCH", "APARCH", "HARCH", "FIGARCH", "MSGARCH"
     ] = "GARCH",
     p: int | list[int] = 1,
     o: int = 0,
@@ -2000,6 +2001,7 @@ def arch_model(
         "harch",
         "constant",
         "egarch",
+        "msgarch",
     )
     known_dist = (
         "normal",
@@ -2060,8 +2062,10 @@ def arch_model(
     elif vol_model == "aparch":
         assert isinstance(p, int)
         v = APARCH(p=p, o=o, q=q)
-    else:  # vol == 'harch'
+    elif vol_model == 'harch':
         v = HARCH(lags=p)
+    else: # vol_model == "msgarch":
+        v = MSGARCH()
 
     if dist_name in ("skewstudent", "skewt"):
         d: Distribution = SkewStudent()