Enhance Hessian handling and improve DSGE model representation

- Added `hessian_strategy` parameter to `run_from_config` and `run_map` functions to allow selection between optimizer and numerical Hessian methods.
- Introduced `build_proposal_from_covariance` function for covariance-based proposal generation.
- Updated `DSGE` class string representation to improve readability and formatting.
- Enhanced LaTeX representation of the DSGE model to ensure clarity in output.
- Added tests to verify the correct LaTeX representation of the DSGE model.

Author: pablo-reyes8

configs/default.yaml

@@ -11,6 +11,7 @@ map:
   method: L-BFGS-B
   hess_step: 0.0001
   tau_scale: 0.5
+  hessian_strategy: auto
   include_jacobian_prior: false
 
 mcmc:

configs/nk_full_yaml.yaml

@@ -131,7 +131,8 @@ map:
   bounds: auto
   include_jacobian_prior: false
   hess_step: 0.0001
-  tau_scale: 0.3
+  tau_scale: 0.25
+  hessian_strategy: auto
 
 mcmc:
   enabled: false

configs/tiny_ar1.yaml

@@ -58,6 +58,7 @@ map:
   include_jacobian_prior: false
   hess_step: 0.0001
   tau_scale: 0.4
+  hessian_strategy: auto
 
 mcmc:
   enabled: false

notebooks/showcase_dsge_octave.ipynb

@@ -130,6 +130,7 @@ def run_from_config(config_path: Path, *, dry_run: bool = False) -> dict[str, An
             "method": map_cfg.get("method", "L-BFGS-B"),
             "hess_step": float(map_cfg.get("hess_step", 1e-4)),
             "tau_scale": float(map_cfg.get("tau_scale", 0.5)),
+            "hessian_strategy": map_cfg.get("hessian_strategy", "auto"),
             "include_jacobian_prior": bool(map_cfg.get("include_jacobian_prior", False)),
         },
         run_mcmc=bool(mcmc_cfg.get("enabled", False)),

scripts/run_pipeline.py

@@ -33,11 +33,15 @@ class ModelSignature:
     n_shocks: int
 
     def __str__(self) -> str:
+        def _label(n: int, singular: str, plural: str) -> str:
+            return f"{n} {singular if n == 1 else plural}"
+
         return (
-            f"{self.n_equations} eqs | "
-            f"{self.n_states} states "
-            f"(+{self.n_leads} leads, {self.n_lags} lags) | "
-            f"{self.n_shocks} shocks")
+            f"{_label(self.n_equations, 'equation', 'equations')} | "
+            f"{_label(self.n_states, 'state', 'states')} "
+            f"({_label(self.n_leads, 'lead', 'leads')}, "
+            f"{_label(self.n_lags, 'lag', 'lags')}) | "
+            f"{_label(self.n_shocks, 'shock', 'shocks')}")
 
 
 def _escape_latex(text: str) -> str:
@@ -174,19 +178,12 @@ def _repr_latex_(self) -> str:
         def _esc(x):
             return _escape_latex(str(x))
 
-        sig = _esc(getattr(self, "signature", ""))
+        sig = str(self.signature).replace("(", "| ").replace(", ", " | ").replace(")", "")
 
         meta = ""
         if getattr(self, "_metadata", None):
             items = [f"{_esc(k)}={_esc(v)}" for k, v in self._metadata.items()]
-            meta = r"\quad\text{\small(" + ", ".join(items) + ")}"
-
-        header = (
-            r"{\normalsize \mathbf{DSGE\ Model}}"
-            r"\\[2pt]"
-            r"{\scriptsize " + sig + r"}"
-            + meta
-            + r"\\[4pt]")
+            meta = r"\\[2pt]{\scriptsize\text{" + ", ".join(items) + r"}}"
 
         rows = []
         for idx, eq in enumerate(self._equations, start=1):
@@ -200,12 +197,16 @@ def _esc(x):
         body = "\n".join(rows)
 
         latex = (
-            r"\begin{aligned}"
-            + header
+            r"\begin{array}{c}"
+            r"{\large\mathbf{DSGE\ Model}}"
+            r"\\[4pt]"
+            r"{\small\text{" + _esc(sig) + r"}}"
+            + meta
+            + r"\\[8pt]"
             + r"\left\{\,\begin{array}{rcll}"
             + body
             + r"\end{array}\right."
-            + r"\end{aligned}")
+            + r"\end{array}")
         return latex
 
 
@@ -303,7 +304,8 @@ def compute(
             allowed_map_keys = {
                 "method",
                 "hess_step",
-                "tau_scale",}
+                "tau_scale",
+                "hessian_strategy",}
 
             filtered_map_kwargs = {k: v for k, v in map_kwargs.items() if k in allowed_map_keys}
 

src/dgse.py

@@ -63,7 +63,7 @@ def numerical_hessian_central(f, x, h=1e-3, ridge=1e-6):
     H = 0.5*(H + H.T) + ridge*np.eye(n)
     return H
 
-def build_proposal_from_hessian(H, tau_scale=0.5, min_eig=1e-8):
+def build_proposal_from_hessian(H, tau_scale=0.5, min_eig=1e-8):
     """
     Toma H ~ ∇²(-logpost)(MAP).
     Devuelve (cov_prop, chol_prop).
@@ -94,7 +94,34 @@ def build_proposal_from_hessian(H, tau_scale=0.5, min_eig=1e-8):
         cov_spd_j = cov_spd + jitter*np.eye(cov_spd.shape[0])
         L = np.linalg.cholesky(cov_spd_j)
 
-    return cov_spd, L
+    return cov_spd, L
+
+
+def build_proposal_from_covariance(cov, tau_scale=0.5, min_eig=1e-8):
+    cov = np.asarray(cov, dtype=float)
+    cov = 0.5 * (cov + cov.T)
+    cov_spd = make_spd_eig(tau_scale * cov, min_eig)
+    try:
+        L = np.linalg.cholesky(cov_spd)
+    except np.linalg.LinAlgError:
+        L = np.linalg.cholesky(cov_spd + 1e-6 * np.eye(cov_spd.shape[0]))
+    return cov_spd, L
+
+
+def _optimizer_inverse_hessian(res, n_params: int):
+    hess_inv = getattr(res, "hess_inv", None)
+    if hess_inv is None:
+        return None
+    try:
+        if hasattr(hess_inv, "todense"):
+            H_inv = np.asarray(hess_inv.todense(), dtype=float)
+        else:
+            H_inv = np.asarray(hess_inv, dtype=float)
+    except Exception:
+        return None
+    if H_inv.shape != (n_params, n_params) or not np.all(np.isfinite(H_inv)):
+        return None
+    return 0.5 * (H_inv + H_inv.T)
 
 
 def run_map(
@@ -114,6 +141,7 @@ def run_map(
     method: str = "L-BFGS-B",
     hess_step: float = 1e-4,
     tau_scale: float = 0.5,
+    hessian_strategy: str = "auto",
     include_jacobian_prior: bool = False , div = 0,
     likelihood_kwargs: Optional[Dict] = None):
 
@@ -147,19 +175,41 @@ def f_obj(theta_work: np.ndarray):
             include_jacobian=include_jacobian_prior,
             likelihood_kwargs=likelihood_kwargs)
 
-    res = minimize(f_obj, theta0, method=method, bounds=bounds)
-
-    theta_map = np.asarray(res.x, dtype=float)
-    H = numerical_hessian_central(f_obj, theta_map, h=hess_step)
-    cov_prop, chol_prop = build_proposal_from_hessian(H, tau_scale=tau_scale, min_eig=1e-8)
+    res = minimize(f_obj, theta0, method=method, bounds=bounds)
+
+    theta_map = np.asarray(res.x, dtype=float)
+    strategy = hessian_strategy.lower()
+    if strategy not in {"auto", "optimizer", "numerical"}:
+        raise ValueError("hessian_strategy debe ser 'auto', 'optimizer' o 'numerical'.")
+
+    H_inv_opt = None
+    if strategy in {"auto", "optimizer"}:
+        H_inv_opt = _optimizer_inverse_hessian(res, theta_map.size)
+
+    if H_inv_opt is not None:
+        cov_prop, chol_prop = build_proposal_from_covariance(
+            H_inv_opt,
+            tau_scale=tau_scale,
+            min_eig=1e-8,
+        )
+        try:
+            H = np.linalg.pinv(H_inv_opt)
+        except np.linalg.LinAlgError:
+            H = np.eye(theta_map.size)
+        hessian_source = "optimizer"
+    else:
+        H = numerical_hessian_central(f_obj, theta_map, h=hess_step)
+        cov_prop, chol_prop = build_proposal_from_hessian(H, tau_scale=tau_scale, min_eig=1e-8)
+        hessian_source = "numerical"
 
 
     return {
         "theta_map": theta_map,
         "neglogpost_map": float(res.fun),
         "success": bool(res.success),
         "message": res.message,
-        "nit": res.nit,
-        "hessian": H,
-        "cov_proposal": cov_prop,
+        "nit": res.nit,
+        "hessian": H,
+        "hessian_source": hessian_source,
+        "cov_proposal": cov_prop,
         "chol_proposal": chol_prop}

src/inference/map.py

@@ -2,3 +2,16 @@ def test_dsge_compatibility_import():
     from src.dsge import DSGE
 
     assert DSGE.__name__ == "DSGE"
+
+
+def test_dsge_latex_repr_has_readable_signature():
+    import sympy as sp
+    from src.dsge import DSGE
+
+    x_t, eps_t = sp.symbols("x_t eps_t")
+    model = DSGE([sp.Eq(x_t, eps_t)], [x_t], eps_t=[eps_t])
+
+    latex = model._repr_latex_()
+
+    assert "1 equation | 1 state | 0 leads | 0 lags | 1 shock" in latex
+    assert "1eqs" not in latex