import json import uuid import numpy as np import zmq from ryp import r, to_py, to_r, options options(to_py_format='numpy') # Load MSGARCH once at import time (equivalent to rpy2's importr('MSGARCH')) r('library(MSGARCH)') class RMSGarchEstimator: """ Wraps an R MSGARCH::CreateSpec object with support for: - arbitrary number of regimes (n_states) - a different variance model / distribution per regime - mixture (do.mix=TRUE) vs Markov-switching (do.mix=FALSE) - ML or MCMC estimation, with MCMC controls exposed fit_and_predict() can return either just the latest-bar summary or the full in-sample regime path so it can be overlaid on the original return series. """ def __init__(self, n_states=2, variance_models="sGARCH", distributions="std", do_mix=False, estimation_method="ML", mcmc_nburn=1000, mcmc_nmesh=2500, mcmc_nthin=1): """ n_states: number of regimes (>= 2) variance_models: single model name (broadcast to all states) or a list of length n_states, e.g. ["sGARCH", "eGARCH", "gjrGARCH"] Any MSGARCH-supported variance model is valid: sGARCH, eGARCH, gjrGARCH, tGARCH, fGARCH, ... distributions: single distribution name (broadcast) or a list of length n_states, e.g. "norm", "std" (Student-t), "ged", "snorm", ... do_mix: True for an iid mixture, False (default) for a Markov-switching model estimation_method: "ML" or "MCMC" mcmc_nburn/mcmc_nmesh/mcmc_nthin: only used when estimation_method == "MCMC" """ if n_states < 2: raise ValueError("n_states must be >= 2") self.n_states = n_states self.variance_models = self._broadcast(variance_models, n_states) self.distributions = self._broadcast(distributions, n_states) self.do_mix = bool(do_mix) self.method = estimation_method self.mcmc_ctr = {"nburn": mcmc_nburn, "nmesh": mcmc_nmesh, "nthin": mcmc_nthin} self._spec_var = f"spec_{uuid.uuid4().hex[:8]}" self._init_model() @staticmethod def _broadcast(val, n): """Allow a single string to be broadcast to all states, or an explicit list.""" if isinstance(val, str): return [val] * n val = list(val) if len(val) != n: raise ValueError(f"Expected {n} entries, got {len(val)}: {val}") return val def _init_model(self): var_vec = "c(" + ", ".join(f'"{v}"' for v in self.variance_models) + ")" dist_vec = "c(" + ", ".join(f'"{d}"' for d in self.distributions) + ")" do_mix_r = "TRUE" if self.do_mix else "FALSE" r(f''' {self._spec_var} <- CreateSpec( variance.spec = list(model = {var_vec}), distribution.spec = list(distribution = {dist_vec}), switch.spec = list(do.mix = {do_mix_r}) ) ''') print(f"[R Engine] Initialized {self.n_states}-state MS-GARCH " f"(variance={self.variance_models}, dist={self.distributions}, " f"do.mix={self.do_mix}, method={self.method})") @staticmethod def _r_matrix_or_none(varname): """Fetch an R object as a numpy array, or None if it's NULL / missing. ryp converts R NULL -> Python None on its own, so this is just a thin, explicitly-typed wrapper around to_py().""" val = to_py(varname) return val @staticmethod def _regime_label(rank, n_states): if n_states == 2: return ["CALM", "CRISIS"][rank] if n_states == 3: return ["CALM", "TRANSITION", "CRISIS"][rank] return f"REGIME_{rank}" # 0 = lowest average volatility ... n_states-1 = highest def fit_and_predict(self, returns, timestamps=None, include_series=False): """ Fits the model on `returns` and returns: - a latest-bar summary (state probabilities, dominant regime, etc.) - if include_series=True, the full in-sample regime path aligned to `returns` / `timestamps`, suitable for overlaying regimes on the original series. Regimes are consistently ordered by ascending average volatility (rank 0 = calmest ... rank n_states-1 = most volatile) rather than by R's arbitrary internal state index, since MSGARCH state labels are not guaranteed to be stable across fits. """ to_r(np.array(returns, dtype=np.float64), 'returns_vec') try: if self.method == "ML": r(f''' fit <- FitML(spec = {self._spec_var}, data = returns_vec) st <- State(fit) predprob <- tryCatch(st$PredProb, error = function(e) NULL) filtprob <- tryCatch(st$FiltProb, error = function(e) NULL) viterbi <- tryCatch(st$Viterbi, error = function(e) NULL) vol <- tryCatch(Volatility(fit), error = function(e) NULL) ''') elif self.method == "MCMC": r(f''' fit <- FitMCMC(spec = {self._spec_var}, data = returns_vec, ctr = list(nburn = {self.mcmc_ctr["nburn"]}, nmesh = {self.mcmc_ctr["nmesh"]}, nthin = {self.mcmc_ctr["nthin"]})) st <- State(fit) predprob <- tryCatch(st$PredProb, error = function(e) NULL) filtprob <- tryCatch(st$FiltProb, error = function(e) NULL) viterbi <- tryCatch(st$Viterbi, error = function(e) NULL) vol <- tryCatch(Volatility(fit), error = function(e) NULL) ''') else: return {"status": "error", "message": f"Unknown estimation_method: {self.method}"} filt_probs = self._r_matrix_or_none('filtprob') pred_probs = self._r_matrix_or_none('predprob') viterbi = self._r_matrix_or_none('viterbi') vol = self._r_matrix_or_none('vol') # Prefer filtered/smoothed probabilities for in-sample regime labeling; # fall back to predicted probabilities if FiltProb isn't available. state_probs = filt_probs if filt_probs is not None else pred_probs if state_probs is None: return {"status": "error", "message": "MSGARCH did not return state probabilities."} n_states = self.n_states if self.method == "MCMC": state_probs = np.mean(state_probs,axis=-2,keepdims=True) pred_probs = np.mean(pred_probs,axis=-2,keepdims=True) viterbi = np.mean(viterbi,axis=-1) viterbi = np.round(viterbi,0) # Order regimes by ascending average volatility across the whole fit, # not just the last bar, so the calm/crisis labeling is stable. if vol is not None and vol.ndim == 2 and vol.shape[1] == n_states: avg_vol = vol.mean(axis=0) rank_order = np.argsort(avg_vol) # rank_order[0] = calmest raw state index else: avg_vol = None rank_order = np.arange(n_states) raw_to_rank = {int(raw): rank for rank, raw in enumerate(rank_order)} latest_probs_ranked = np.reshape(state_probs[-1],state_probs.shape[-1]) dominant_state_raw = int(np.argmax(latest_probs_ranked)) dominant_rank = raw_to_rank[dominant_state_raw] summary = { "status": "success", "method": self.method, "n_states": n_states, "variance_models": self.variance_models, "distributions": self.distributions, "do_mix": self.do_mix, "state_probabilities": latest_probs_ranked.tolist(), # ordered calm -> turbulent "dominant_regime": dominant_rank, "regime_label": self._regime_label(dominant_rank, n_states), "state_volatility": (avg_vol[rank_order].tolist() if avg_vol is not None else None), } # Backward-compatible fields for existing 2-state consumers (e.g. CSignalGarchTouch.mqh) if n_states == 2: summary["p_calm"] = float(latest_probs_ranked[0]) summary["p_crisis"] = float(latest_probs_ranked[1]) summary["current_state"] = "CRISIS" if latest_probs_ranked[1] > 0.5 else "CALM" summary["vol_calm"] = summary["state_volatility"][0] if avg_vol is not None else None summary["vol_crisis"] = summary["state_volatility"][1] if avg_vol is not None else None if include_series: T = state_probs.shape[0] ranked_probs_series = state_probs.reshape((state_probs.shape[0],state_probs.shape[-1])) dominant_per_bar_raw = np.argmax(state_probs, axis=2) dominant_rank_per_bar = np.array([raw_to_rank[int(r_[0])] for r_ in dominant_per_bar_raw]) vol_series_ranked = None if vol is not None and vol.ndim == 2 and vol.shape[1] == n_states: vol_series_ranked = vol[:, rank_order] viterbi_rank_series = None if viterbi is not None: # MSGARCH's Viterbi path is 1-indexed R state labels viterbi_raw = np.array(viterbi).astype(int).flatten() - 1 viterbi_rank_series = np.array([raw_to_rank[int(r_)] for r_ in viterbi_raw]) returns_arr = np.array(returns, dtype=np.float64) series = [] for i in range(T): entry = { "index": i, "timestamp": timestamps[i] if timestamps is not None and i < len(timestamps) else None, "return": float(returns_arr[i]) if i < len(returns_arr) else None, "state_probs": ranked_probs_series[i].tolist(), "regime": int(dominant_rank_per_bar[i]), "regime_label": self._regime_label(int(dominant_rank_per_bar[i]), n_states), } if vol_series_ranked is not None: entry["state_vol"] = vol_series_ranked[i].tolist() if viterbi_rank_series is not None and i < len(viterbi_rank_series): entry["viterbi_regime"] = int(viterbi_rank_series[i]) series.append(entry) summary["series"] = series return summary except Exception as e: return {"status": "error", "message": str(e)} DEFAULT_CONFIG_KEYS = ( "n_states", "variance_models", "distributions", "do_mix", "estimation_method", "mcmc_nburn", "mcmc_nmesh", "mcmc_nthin", ) def start_server(port=5555, default_config=None): """ ZMQ REP server. Two request types (JSON): 1) {"action": "configure", ...any of DEFAULT_CONFIG_KEYS...} Updates the server's default MSGARCH configuration for future "predict" calls that don't supply their own per-call "config". 2) {"action": "predict", "returns": [...], # required, >= 50 observations "timestamps": [...], # optional, same length as returns "include_series": true/false, # optional, default false — set true # to get the full regime-overlay series "config": {...}} # optional one-off override of any # DEFAULT_CONFIG_KEYS, does not persist Estimators are cached by their full config so repeated configs (the common case) don't pay the R CreateSpec cost every request. """ default_config = default_config or {} context = zmq.Context() socket = context.socket(zmq.REP) socket.bind(f"tcp://127.0.0.1:{port}") estimators_cache = {} def get_estimator(cfg): key = json.dumps(cfg, sort_keys=True) if key not in estimators_cache: estimators_cache[key] = RMSGarchEstimator(**cfg) return estimators_cache[key] current_default_cfg = { "n_states": default_config.get("n_states", 2), "variance_models": default_config.get("variance_models", "sGARCH"), "distributions": default_config.get("distributions", "std"), "do_mix": default_config.get("do_mix", False), "estimation_method": default_config.get("estimation_method", "ML"), "mcmc_nburn": default_config.get("mcmc_nburn", 1000), "mcmc_nmesh": default_config.get("mcmc_nmesh", 2500), "mcmc_nthin": default_config.get("mcmc_nthin", 1), } get_estimator(current_default_cfg) # warm the default spec on startup print(f"[ZMQ] MS-GARCH Server active on port {port} | default config: {current_default_cfg}") while True: try: msg_str = socket.recv_string() data = json.loads(msg_str) action = data.get("action") if action == "configure": for k in DEFAULT_CONFIG_KEYS: if k in data: current_default_cfg[k] = data[k] try: get_estimator(current_default_cfg) # validate + warm resp = {"status": "success", "message": "Default configuration updated.", "config": current_default_cfg} except Exception as e: resp = {"status": "error", "message": f"Invalid configuration: {e}"} elif action == "predict": returns = data.get("returns", []) timestamps = data.get("timestamps") include_series = bool(data.get("include_series", False)) cfg_override = data.get("config") or {} cfg = dict(current_default_cfg) cfg.update({k: v for k, v in cfg_override.items() if k in DEFAULT_CONFIG_KEYS}) if len(returns) < 50: resp = {"status": "error", "message": "Min 50 return observations required for R fit."} else: try: estimator = get_estimator(cfg) resp = estimator.fit_and_predict(returns, timestamps=timestamps, include_series=include_series) except Exception as e: resp = {"status": "error", "message": f"Invalid configuration: {e}"} else: resp = {"status": "error", "message": "Invalid action. Use 'predict' or 'configure'."} socket.send_string(json.dumps(resp)) except Exception as e: socket.send_string(json.dumps({"status": "error", "message": str(e)})) if __name__ == "__main__": # Example: start with a 3-state model, GJR-GARCH in every regime, Student-t innovations start_server( port=5555, default_config={ "n_states": 2, "variance_models": "sGARCH", "distributions": "std", "do_mix": False, "estimation_method": "ML", }, )