Skip to content

🇬🇧 English · 🇪🇸 Español

Cuaderno de bocetos DSGE: modelos desde ecuaciones, no matrices

puremacro.dsge.klein_solve siempre ha resuelto modelos lineales de expectativas racionales sin necesidad de Dynare ni de un compilador externo — pero requiere las matrices A, B y C, y obtenerlas implica diferenciar las condiciones de equilibrio a mano. Esa derivación algebraica es precisamente el paso donde los investigadores cometen errores y el paso en el que una tableta ofrece menos facilidades: sin sistema algebraico computacional, sin MATLAB y sin margen para el error.

dsge.build elimina esa barrera. Escriba las condiciones de equilibrio tal como aparecen en el artículo:

import numpy as np
from puremacro import dsge

def eqs(xp, x, e, p):
    # xp = t+1 (adelantos), x = t (actuales), e = choques, p = parámetros
    return [
        1 / x.c - p.beta * (p.alpha * xp.z * xp.k ** (p.alpha - 1)) / xp.c,
        x.c + xp.k - x.z * x.k ** p.alpha,
        xp.z - x.z ** p.rho * np.exp(e.eps),
    ]

m = dsge.build(
    eqs,
    variables=["c", "k", "z"],
    states=["k", "z"],
    shocks=["eps"],
    params=dict(alpha=0.33, beta=0.98, rho=0.9),
    guess=dict(c=0.5, k=0.1, z=1.0),
)
print(m.summary())
m.irf("eps", horizon=20)     # DataFrame: horizontes × variables

xp, x, e y p admiten acceso por atributos (x.k), indexación por cadenas (x["k"]), indexación posicional y desempaquetado de tuplas. De este modo, el modelo se lee en el código exactamente igual que en el papel. Un nombre mal escrito genera inmediatamente un error informativo listando las variables declaradas.

Qué devuelve el resolvedor

LinearModel es una dataclass congelada que almacena el estado estacionario, las unidades, las matrices en forma de Klein y la KleinSolution subyacente:

Miembro Significado
.steady_state pandas.Series, verificado contra las ecuaciones en lugar de asumido
.units "log" o "level" por cada variable
.irf(shock, horizon, size) Funciones de respuesta al impulso (horizontes × variables)
.simulate(periods, sigma, seed) Simulación estocástica lineal
.policy() Reglas de decisión en una tabla etiquetada
.decision_rules() / .dynare_dr Objeto DynareDR idéntico a oo_.dr de Dynare (ghx, ghu, ys)
.theoretical_moments() Momentos analíticos, correlaciones cruzadas, autocorrelaciones y FEVD (stoch_simul)
.fevd(horizons) Descomposición analítica de la varianza del error de pronóstico
.solution La KleinSolution resuelta — G, F, N, L
.A, .B, .C Las matrices producidas por la diferenciación de paso complejo

Los estados estacionarios se calculan numéricamente a partir de guess=, o bien se suministran mediante steady_state= y son verificados rigurosamente contra las condiciones de equilibrio (max|f(ss, ss, 0)| <= tol). La log-linealización opera variable por variable, recurriendo a desviaciones en niveles cuando el estado estacionario no es estrictamente positivo.

Diferenciación por paso complejo y sus restricciones

Los jacobianos se calculan mediante \(\operatorname{Im} f(x + ih) / h\) con \(h = 10^{-20}\) — derivadas con precisión de máquina en una única evaluación por argumento, sin compensaciones de tamaño de paso y sin errores de cancelación numérica.

La única restricción matemática es que este método es exacto si y sólo si la función residual es analítica, y falla silenciosamente si no lo es: \(\operatorname{Im}\) evaluada a través de abs() es idénticamente cero, devolviendo derivadas nulas sin generar excepciones. Cuatro construcciones rompen la analiticidad: abs, min/max, comparaciones condicionales que ramifican según valores perturbados y conversiones explícitas a float() o np.real() que descartan la parte imaginaria.

Por ello, build verifica de forma cruzada cada jacobiano de paso complejo contra diferencias finitas en una dirección aleatoria, generando un ModelError si discrepan. Para modelos con restricciones ocasionalmente activas, pase method="central" con precisión de \(\approx 10^{-8}\).

Temporización y orden temporal

Las ecuaciones se plantean como \(\mathbb{E}_t f(z_{t+1}, z_t, u_t) = 0\), respetando la convención de klein_solve. Un proceso exógeno estándar \(z_{t+1} = \rho z_t + \varepsilon_{t+1}\) traslada el estado al período siguiente, por lo que todo estado es cero en la fila \(h=0\) de una FRI e impacta en \(h=1\).

Por el contrario, los controles con visión de futuro (forward-looking) habitualmente sí reaccionan en \(h=0\), ya que la innovación entra en el conjunto de información de los agentes y estos descuentan el valor futuro en \(t+1\). Ese salto contemporáneo corresponde a la matriz de carga \(L\) de Klein.

Condiciones de Blanchard-Kahn

La determinación no es una premisa impuesta sino una propiedad comprobada durante la resolución. Si se viola el principio de Taylor en un modelo nuevo keynesiano, el resolvedor levanta una excepción explícita:

solve(phi_pi=0.9)
# BlanchardKahnError: Blanchard-Kahn indeterminacy: 2 unstable generalised
# eigenvalues vs 3 forward-looking variables

Pase strict=False si prefiere el comportamiento permisivo de Sims/gensys (matrices nulas con bandera de error eu).

Validación analítica

El modelo neoclásico de crecimiento con depreciación completa y utilidad logarítmica admite solución analítica exacta: \(k' = \alpha\beta z k^\alpha\), \(c = (1-\alpha\beta) z k^\alpha\), con \(G = [[\alpha, 1], [0, \rho]]\) y \(F = [\alpha, 1]\). build reproduce cada celda de \(G\), \(F\), \(N\) y \(L\) con un error inferior a \(10^{-9}\) y el estado estacionario hasta \(10^{-12}\).

Compatibilidad con Dynare: Reglas de Decisión y Momentos

Para reproducir fielmente la salida de stoch_simul(order=1) de Dynare, LinearModel proporciona:

1. Reglas de decisión (m.decision_rules() / m.dynare_dr)

Mapea la solución de Klein a la representación canónica de Dynare: $\(y_t = y^* + g_x (x_{t-1} - x^*) + g_u u_t\)$

dr = m.decision_rules()
print(dr.summary())
Imprime la tabla exacta POLICY AND TRANSITION FUNCTIONS con la constante, los estados rezagados (k(-1), z(-1)) y los choques estructurales.

2. Momentos teóricos analíticos (m.theoretical_moments())

Resuelve la ecuación discreta de Lyapunov \(\Sigma_x = G \Sigma_x G' + N \Sigma_u N'\) para obtener los momentos no condicionados analíticos exactos sin recurrir a simulaciones de Monte Carlo:

mom = m.theoretical_moments(lags=5)
print(mom.summary())
Genera los 4 bloques estándar de Dynare: - THEORETICAL MOMENTS: Media, desviación estándar y varianza de cada variable. - MATRIX OF CORRELATIONS: Matriz completa de correlaciones cruzadas. - COEFFICIENTS OF AUTOCORRELATION: Autocorrelaciones teóricas exactas para los retardos 1 a 5. - VARIANCE DECOMPOSITION: Participaciones porcentuales de la descomposición de varianza del error de pronóstico (FEVD) en horizontes finitos (1, 4, 8, 16, 32) y asintótico.

3. Interfaz canónica de adelantos y retardos (dsge.build_dynare)

Permite especificar modelos en la representación dinámica canónica de Dynare: $\(\mathbb{E}_t [ f(y_{t+1}, y_t, y_{t-1}, u_t; \theta) ] = 0\)$

from puremacro import dsge

def rbc(lead, curr, lag, shocks, p):
    return [
        curr.c**(-p.gamma) - p.beta * lead.c**(-p.gamma) * (p.alpha * np.exp(lead.a) * curr.k**(p.alpha - 1) + 1 - p.delta),
        curr.k - (np.exp(curr.a) * lag.k**p.alpha - curr.c + (1 - p.delta) * lag.k),
        curr.a - (p.rho * lag.a + shocks.eps),
    ]

m = dsge.build_dynare(
    rbc,
    variables=["k", "a", "c"],
    shocks=["eps"],
    params=dict(alpha=0.3, beta=0.99, delta=0.025, gamma=1.0, rho=0.8),
    guess=dict(k=38.0, a=0.0, c=2.0),
)
# ¡Las variables de estado ('k', 'a') y control ('c') se clasifican automáticamente!

4. Analizador nativo de archivos .mod de Dynare (dsge.load_mod)

Ejecute archivos .mod directamente en Python puro sin dependencias de MATLAB u Octave:

from puremacro import dsge

mod_text = """
var c k a;
varexo eps;
parameters alpha beta delta gamma rho;

alpha = 0.30;
beta  = 0.99;
delta = 0.025;
gamma = 1.0;
rho   = 0.80;

model;
  c^(-gamma) = beta * c(+1)^(-gamma) * (alpha * exp(a(+1)) * k^(alpha - 1.0) + 1.0 - delta);
  k = exp(a) * k(-1)^alpha - c + (1.0 - delta) * k(-1);
  a = rho * a(-1) + eps;
end;

initval;
  k = 38.0;
  a = 0.0;
  c = 2.0;
end;
"""

# Cargar desde texto o desde una ruta de archivo
m = dsge.load_mod(mod_text)
print(m.decision_rules().summary())
print(m.theoretical_moments().summary())

4b. Nivel 0: El Front End — Preprocesador de Macros, Parser AST y Diferenciación Simbólica (v2.7.0)

A partir de puremacro 2.7.0, el motor DSGE incorpora un pipeline frontal de nivel de producción diseñado para compatibilidad nativa con modelos de Dynare publicados y sin ninguna dependencia externa:

Directivas de Macro Dynare

puremacro.dsge preprocesa de forma nativa las directivas macro de Dynare antes de la tokenización: - @#define VAR = EXPR: define constantes enteras, reales, booleanas, cadenas o arreglos. - @#for VAR in EXPR ... @#endfor: iteración sobre rangos numéricos (1:N) o listas literales (["c", "k"]). - @#if EXPR ... @#elif EXPR ... @#else ... @#endif: compilación condicional según parámetros o banderas de régimen. - @#include "ruta/al/archivo.mod": inclusión recursiva de archivos con resolución de rutas relativas. - @{EXPR}: interpolación de expresiones y valores directamente en ecuaciones o parámetros.

mod_macro = """
@#define CALIBRATION = "baseline"
@#define SECTORS = ["agr", "mfg", "srv"]

var
@#for s in SECTORS
  y_@{s} c_@{s}
@#endfor
;
...
"""
m = dsge.load_mod(mod_macro)

Árbol de Sintaxis Abstracta (AST) y Parser Descendente Recursivo

  • Grafo Acíclico Dirigido Inmutable (DAG): Las ecuaciones se transforman en un grafo inmutable de expresiones (Const, Var, Param, UnaryOp, BinOp, Call), eliminando manipulaciones frágiles de texto y errores por colisiones de subcadenas.
  • Variables Locales #: Admite declaraciones como #MU = c^(-gamma); y #cbetabar = cbeta*cgamma^(-csigma);. Como en Dynare (Manual de Referencia, §4.5), una variable local es una abreviatura simbólica: su expresión se sustituye allí donde se usa, se desplaza en el tiempo cuando aparece como MU(+1) o MU(-1), y nunca se evalúa a un número al cargar el archivo. Por eso una local construida a partir de parámetros los sigue cada vez que el modelo se vuelve a resolver: load_mod(..., params={...}), osr, los widgets, los gradientes analíticos, SMC y model.estimate() ven todos el valor nuevo (el constebeta de SW07 entra en el modelo solo a través de variables locales). dag.evaluate_locals(params) devuelve sus valores numéricos. Hasta la 4.3.0 inclusive, una local que dependía solo de parámetros quedaba congelada en la calibración del archivo, de modo que cada nueva resolución conservaba su valor de carga.
  • Bloques que dependen de parámetros al sobrescribirlos: load_mod(..., params={...}) vuelve a evaluar steady_state_model e initval con los parámetros combinados (Dynare vuelve a ejecutar steady_state_model en cada vector de parámetros). Las asignaciones de primer nivel como beta = 1/(1+r); se evalúan una sola vez al leer el archivo, como en Dynare, y no cambian si se sobrescribe r.
  • Cualquier nombre de Dynare: un parámetro o una variable con nombre de palabra reservada de Python (lambda, yield) se analiza y se resuelve.
  • Operadores Especiales: Soporte para STEADY_STATE(x) y diff(x). EXPECTATION(k)(x) se analiza sintácticamente, pero su residuo compilado lanza NameError, así que un modelo que lo use todavía no puede resolverse.
  • Funciones Trascendentes: Análisis sintáctico y diferenciación analítica exacta para exp, log, sin, cos, tan, normcdf, normpdf, erf, abs y sign.

Diferenciación Simbólica y Eliminación de Subexpresiones Comunes (CSE)

  • Jacobianos de primer orden (\(A_+, A_0, A_-, B_u\)) y matriz hessiana dinámica de segundo orden (\(H_f\)) diferenciados analíticamente con cero error de aproximación numérica.
  • La Eliminación de Subexpresiones Comunes (CSE) optimiza el grafo eliminando entre un 30% y un 70% de operaciones algebraicas redundantes, compilando derivados dinámicos en evaluadores de alto rendimiento.

Solucionador Rápido Schur Sylvester

  • La ecuación de curvatura de segundo orden \((A_0 + A_+ g_x P_s) g_{xx} + A_+ g_{xx} (h_x \otimes h_x) = -K_{xx}\) se resuelve mediante la descomposición compleja de Schur \(h_x = U_h T_h U_h^H\) con almacenamiento en caché de factorizaciones LU diagonales y propagación hacia adelante.
  • Resuelve modelos a gran escala como Smets-Wouters (2007) en 0.028 segundos (frente a 4.6s anteriores), reduciendo el consumo de memoria de 648 MB a menos de 2 MB.

Introspección del Modelo y Exportación a LaTeX

  • Diagnóstico Estructural: m.model_info() o dsge.model_info(dag) proporciona un reporte detallado de incidencias temporales, clasificación de variables y dispersión del jacobiano.
  • Ecuaciones en LaTeX para Publicación: m.write_latex_dynamic_model("modelo.tex") genera bloques \begin{align}...\end{align} listos para su inclusión en manuscritos académicos.

5. Perturbación de segundo orden con poda (Pruning)

Resuelva aproximaciones cuadráticas estables siguiendo el algoritmo de Schmitt-Grohé y Uribe (2004) y Kim, Kim, Schaumburg y Sims (2008):

# Segundo orden: pase order=2 a load_mod / build_dynare, o vuelva a resolver un modelo existente
sol_2nd = m.solve(order=2)              # PrunedDSGESolution (Kim-Kim-Schaumburg-Sims pruning)
print(sol_2nd.oo_dr.summary())          # ghx, ghu, ghxx, ghxu, ghuu, ghs2 en formato Dynare

Descompone el espacio de estados en componentes de primer y segundo orden: $\(x_t^{(1)} = G x_{t-1}^{(1)} + N u_t\)$ $\(x_t^{(2)} = G x_{t-1}^{(2)} + \frac{1}{2} H_{xx} (x_{t-1}^{(1)} \otimes x_{t-1}^{(1)}) + \frac{1}{2} H_{\sigma\sigma} \sigma^2\)$ $\(y_t^{(1)} = F x_t^{(1)} + L u_t\)$ $\(y_t^{(2)} = F x_t^{(2)} + \frac{1}{2} G_{xx} (x_t^{(1)} \otimes x_t^{(1)}) + \frac{1}{2} G_{\sigma\sigma} \sigma^2\)$

Capacidades destacadas: - Simulación incondicionalmente estable: sol_2nd.simulate(periods=200) evita las trayectorias explosivas características de las aproximaciones cuadráticas sin poda. Su sigma multiplica la escala de los choques ya declarada en shock_cov, de modo que sigma=1 (el valor por defecto) es la calibración. - Respuestas al impulso generalizadas (GIRF): sol_2nd.girf(shock="eps", size=0.01, horizon=20) evalúa impactos dependientes del estado. - Media ergódica y estado estacionario con riesgo: sol_2nd.ergodic_mean() (alias stochastic_steady_state()) es la media incondicional de la solución podada; sol_2nd.risky_steady_state() es el punto fijo sin choques, donde solo actúa el término de riesgo ½ g_σσ σ²; sol_2nd.risk_decomposition() descompone la media en ese término de riesgo y los términos de curvatura ½ g_xx vec(Ω) y ½ g_uu vec(σ²Σ_u) (véase Perturbación de orden superior, §2.4). Ambos pueden tener signos opuestos: en el modelo RBC de esa sección la media ergódica del capital está por encima de su estado estacionario determinista y el estado estacionario con riesgo, por debajo.

6. Herramienta de línea de comandos puremacro-dynare CLI

puremacro 2.2.0 incluye el comando de terminal puremacro-dynare (puremacro.dsge.cli), permitiendo resolver modelos desde scripts de consola o terminales:

# Resolución básica y visualización de reglas de política
puremacro-dynare rbc.mod

# Resolver a 2do orden con poda y generar tabla FEVD
puremacro-dynare rbc.mod --order 2 --fevd

# Exportar tablas listas para publicación (Markdown, LaTeX, Typst) y guardar gráficos FRI
puremacro-dynare rbc.mod --irf 20 --format all --plot

# Descomposición histórica de choques a partir de un archivo CSV de datos observados
puremacro-dynare rbc.mod --shock-decomp datos_macro.csv --plot

7. OccBin: Restricciones Ocasionalmente Activas (ZLB)

puremacro.dsge.occbin implementa el algoritmo lineal por tramos de Guerrieri e Iacoviello (2015) para modelos con cotas inferiores de tasa de interés cero (Zero Lower Bound, ZLB) o restricciones de colateral.

Calcula la trayectoria de alternancia de regímenes entre el régimen de referencia \(M_1\) y el régimen restringido \(M_2\) mediante recursión hacia atrás:

import numpy as np
from puremacro.dsge import build_dynare, solve_occbin, OccBinConstraint

# Modelo neokeynesiano de tres ecuaciones. El régimen restringido fija la tasa
# nominal en el límite inferior cero (r = -r_ss en desviaciones del estado estacionario).
params = {"beta": 0.99, "sigma": 1.0, "kappa": 0.1, "phi_pi": 1.5, "phi_y": 0.125, "rho_g": 0.8, "r_ss": 0.01}
variables = ["y", "pi", "r", "g"]
shocks = ["eps_r", "eps_g"]
steady_state = {v: 0.0 for v in variables}

def nk_taylor(lead, curr, lag, shocks_v, p):
    return [
        curr.y - lead.y + (curr.r - lead.pi) / p.sigma - curr.g,          # dynamic IS
        curr.pi - p.beta * lead.pi - p.kappa * curr.y,                     # NK Phillips curve
        curr.r - p.phi_pi * curr.pi - p.phi_y * curr.y - shocks_v.eps_r,   # Taylor rule
        curr.g - p.rho_g * lag.g - shocks_v.eps_g,                         # demand shock process
    ]

def nk_zlb(lead, curr, lag, shocks_v, p):
    return [
        curr.y - lead.y + (curr.r - lead.pi) / p.sigma - curr.g,
        curr.pi - p.beta * lead.pi - p.kappa * curr.y,
        curr.r - (-p.r_ss),                                                # rate pegged at the ZLB
        curr.g - p.rho_g * lag.g - shocks_v.eps_g,
    ]

m_ref = build_dynare(nk_taylor, variables=variables, shocks=shocks, params=params, steady_state=steady_state)
# Un régimen con tasa fija es indeterminado por sí solo; OccBin solo necesita sus
# jacobianos, por eso el régimen alternativo se construye con strict=False.
m_zlb = build_dynare(nk_zlb, variables=variables, shocks=shocks, params=params,
                     steady_state=steady_state, check_steady_state=False, strict=False)

# La restricción se activa cuando la tasa sombra cae por debajo del límite
constraint = OccBinConstraint(variable="r", threshold=-params["r_ss"], operator="<")

horizon = 30
shock_seq = np.zeros((horizon, len(shocks)))
shock_seq[0, shocks.index("eps_g")] = -0.04        # choque de demanda contractivo en t=0

res_occbin = solve_occbin(m_ref, m_zlb, constraint, shock_sequence=shock_seq, horizon=horizon)
print(res_occbin.summary())
print("Períodos en el límite inferior cero:", res_occbin.binding_periods)
print("Régimen por período (1 = restringido):", res_occbin.regimes)
res_occbin.plot()

8. Simulación no lineal y previsión perfecta

puremacro.dsge.perfect_foresight implementa el método de relajación apilada de Newton-Raphson de Boucekkine (1995) y Juillard (1996) para resolver transiciones deterministas no lineales exactas:

import numpy as np
from puremacro.dsge import solve_perfect_foresight

# Modelo de Ramsey determinista en niveles:
#   1/c_t = beta / c_{t+1} * (alpha A_t k_t^(alpha-1) + 1 - delta)
#   k_t   = A_t k_{t-1}^alpha + (1 - delta) k_{t-1} - c_t
alpha, beta, delta = 0.33, 0.96, 0.10
r_ss = 1.0 / beta - (1.0 - delta)
k_ss = (alpha / r_ss) ** (1.0 / (1.0 - alpha))
c_ss = k_ss ** alpha - delta * k_ss

def ramsey(y_plus, y_curr, y_lag, exo):
    c_p, k_p = y_plus
    c, k = y_curr
    c_m, k_m = y_lag
    A = float(np.ravel(exo)[0])
    return [
        1.0 / c - beta / c_p * (alpha * A * k ** (alpha - 1.0) + 1.0 - delta),
        k - (A * k_m ** alpha + (1.0 - delta) * k_m - c),
    ]

# Un alza de +5% en la PTF en el período 5, anunciada en t=0: el consumo salta antes de que llegue el choque
n_periods = 100
tfp_path = np.ones(n_periods)
tfp_path[4] = 1.05

pf_res = solve_perfect_foresight(
    ramsey,
    y_init=np.array([c_ss, k_ss]),
    y_ss=np.array([c_ss, k_ss]),
    exogenous_path=tfp_path,
    n_periods=n_periods,
    variable_names=["c", "k"],
)
print(pf_res.summary())
pf_res.plot()

9. Estimación Bayesiana de DSGE por MCMC

puremacro.dsge.bayesian proporciona un flujo bayesiano completo sin dependencias compiladas externas: 1. Distribuciones a priori (BetaPrior, GammaPrior, InvGammaPrior, NormalPrior, UniformPrior). 2. Búsqueda de moda mediante L-BFGS-B o Nelder-Mead. 3. Inversión del hessiano numérico en la moda para la covarianza de propuesta de Laplace. 4. Muestreador adaptativo Random-Walk Metropolis-Hastings (RWMH). 5. Diagnósticos de convergencia \(\hat{R}\) dividida (Gelman-Rubin) y prueba espectral de Geweke, con la densidad espectral en cero estimada a partir de un modelo AR ajustado, como en coda de R.

import numpy as np
from puremacro.dsge import estimate_dsge_bayesian
from puremacro.dsge.priors import BetaPrior, InvGammaPrior
from puremacro.state_space import StateSpaceModel, kalman_filter

# Observable: y_t = rho y_{t-1} + sigma eps_t, T = 300 (sustituto de la verosimilitud de su modelo)
rng = np.random.default_rng(42)
y = np.zeros(300)
for t in range(1, 300):
    y[t] = 0.7 * y[t - 1] + 0.4 * rng.standard_normal()
y = y[:, None]

def log_likelihood(params):
    rho, sigma = (params["rho"], params["sigma"]) if isinstance(params, dict) else (params[0], params[1])
    ssm = StateSpaceModel(T=np.array([[rho]]), Z=np.array([[1.0]]), R=np.array([[1.0]]),
                          Q=np.array([[sigma ** 2]]), H=np.array([[1e-6]]))
    return kalman_filter(y, ssm)["loglik"]

priors = {
    "rho": BetaPrior(mean=0.6, std=0.15, lb=0.01, ub=0.99),
    "sigma": InvGammaPrior(mean=0.3, std=2.0, lb=0.01, ub=3.0),
}

bayes_res = estimate_dsge_bayesian(
    log_likelihood, priors,
    initial_params=np.array([0.6, 0.3]),
    n_draws=1000, n_burn=200, n_chains=2, seed=42,
)
print(bayes_res.summary())

# Densidades a priori y a posteriori, y una tabla lista para publicación
bayes_res.plot_priors_posteriors()
print(bayes_res.to_latex())

10. Descomposición Histórica de Choques y FEVD

# 1. FEVD analítica en los horizontes [1, 4, 8, 16, 40]
fevd_res = m.fevd_result(horizons=[1, 4, 8, 16, 40])
print(fevd_res.summary())
print(fevd_res.to_latex())

# 2. Descomposición histórica de choques de una muestra (aquí simulada del modelo;
#    pase sus datos observados con una columna por variable)
data = m.simulate(periods=80, seed=1)
decomp_res = m.shock_decomposition(data)
print(decomp_res.summary())

# Gráfico de barras apiladas de las contribuciones históricas de los choques
decomp_res.plot(variable="c")