API Reference: Simulation & IRFs¶

def irfs(
    model: AbstractModel,
    dr: RecursiveDecisionRule | PerturbationSolution,
    type: IRFType = "log-deviation",
    T: int = 40,
) -> IRFSimulation:
    """Impulse response function simulation in response to shocks on each exogenous variable."""
    exo_names = list(model.symbols["exogenous"])
    endo_names = list(dr.symbols["endogenous"])
    n_endo = len(endo_names)

    if len(exo_names) == 0:
        data = np.zeros((0, T + 1, n_endo), dtype=float)
    else:
        arrays = [_compute_irf_array(dr, i, T=T) for i in range(len(exo_names))]
        data = np.stack(arrays, axis=0)

    x0 = dr.x0 if dr.x0 is not None else np.zeros(n_endo, dtype=float)
    return IRFSimulation(
        data=data,
        variables=endo_names,
        shocks=exo_names,
        steady_state=x0,
        units=type,
        canonical_units="deviation",
    )

def simulate(
    dr: RecursiveDecisionRule | PerturbationSolution,
    T: int = 40,
    shocks: dict[str, dict[int, float]] | np.ndarray | None = None,
    initial_states: dict[str, float] | np.ndarray | None = None,
    N: int = 1,
    mode: SimulateMode = "random",
    units: UnitsType = "deviation",
    rng: np.random.Generator | int | None = None,
) -> SimulationResult:
    """Simulate the evolution of endogenous variables from a solved decision rule.

    Parameters
    ----------
    dr : RecursiveDecisionRule | PerturbationSolution
        Solved linearized model.
    T : int, default 40
        Time horizon over which the simulation is performed.
    shocks : dict[str, dict[int, float]] | np.ndarray | None, optional
        Preset shocks to force at specific dates.
    initial_states : dict[str, float] | np.ndarray | None, optional
        Preset initial state at date 0.
    N : int, default 1
        Number of random trajectories to simulate when ``mode='random'``.
    mode : {'auto', 'irf', 'random', 'transition', 'deterministic'}, default 'random'
        Simulation mode.
    units : {'deviation', 'level', 'percent', 'log-deviation'}, default 'deviation'
        Default output units for the returned ``SimulationResult``.
    rng : numpy.random.Generator | int | None, optional
        Random number generator, or a seed, used to draw the shocks when
        ``mode='random'``. ``None`` draws from a fresh unseeded generator.

    Returns
    -------
    SimulationResult
        ``IRFSimulation`` when ``mode='irf'``, or ``RandomSimulation`` of shape ``(N, T+1, V)``.
    """
    generator = np.random.default_rng(rng)
    if mode == "irf":
        model = getattr(dr, "_model", None) or getattr(
            getattr(dr, "decision_rule", None), "_model", None
        )
        if model is None:
            raise ValueError("Cannot compute IRFs without an attached model.")
        return irfs(model, dr, type=units, T=T)

    X = dr.X
    Y = dr.Y
    Σ = dr.Σ
    assert X.shape is not None
    assert Y.shape is not None

    n_endo = X.shape[1]
    n_exo = Y.shape[1]
    m0 = np.zeros(n_exo)

    endo_names = list(dr.symbols["endogenous"])
    exo_names = list(dr.symbols["exogenous"])

    model = getattr(dr, "_model", None) or getattr(
        getattr(dr, "decision_rule", None), "_model", None
    )
    context = getattr(model, "context", {}) if model is not None else {}
    model_values = context.get("values", {})
    steady_states = context.get("steady_states", {})

    # 1. Initial shock at t=0
    eps_0 = np.zeros(n_exo)
    for j, name in enumerate(exo_names):
        if (
            shocks is not None
            and isinstance(shocks, dict)
            and name in shocks
            and 0 in shocks[name]
        ):
            val = shocks[name][0]
            ss_val = steady_states.get(name, 0.0)
            eps_0[j] = val - ss_val if not np.isnan(ss_val) else val
        elif name in model_values and 0 in model_values[name]:
            val = model_values[name][0]
            ss_val = steady_states.get(name, 0.0)
            eps_0[j] = val - ss_val if not np.isnan(ss_val) else val

    # 2. Initial state at t=0
    v0 = Y @ eps_0

    if initial_states is not None:
        if isinstance(initial_states, dict):
            for i, name in enumerate(endo_names):
                if name in initial_states:
                    val = initial_states[name]
                    ss_val = steady_states.get(name, 0.0)
                    v0[i] = val - ss_val if not np.isnan(ss_val) else val
        elif isinstance(initial_states, (np.ndarray, list)):
            v0 = np.array(initial_states, dtype=float).copy()
    else:
        for i, name in enumerate(endo_names):
            if name in model_values and 0 in model_values[name]:
                val = model_values[name][0]
                ss_val = steady_states.get(name, 0.0)
                v0[i] = val - ss_val if not np.isnan(ss_val) else val

    n_draws = max(1, int(N))
    trajectories = np.zeros((n_draws, T + 1, n_endo), dtype=float)
    trajectories[:, 0, :] = v0[None, :]

    # 3. Dynamic simulation for t = 1 to T across N draws
    for draw_idx in range(n_draws):
        curr = v0.copy()
        for t in range(1, T + 1):
            if shocks is not None and isinstance(shocks, np.ndarray):
                if shocks.ndim == 3:
                    e = shocks[draw_idx, t - 1, :].copy()
                else:
                    e = shocks[t - 1, :].copy()
            else:
                if n_exo > 0:
                    e = generator.multivariate_normal(m0, Σ)
                else:
                    e = np.zeros(0, dtype=float)
                for j, name in enumerate(exo_names):
                    if (
                        shocks is not None
                        and isinstance(shocks, dict)
                        and name in shocks
                        and t in shocks[name]
                    ):
                        val = shocks[name][t]
                        ss_val = steady_states.get(name, 0.0)
                        e[j] = val - ss_val if not np.isnan(ss_val) else val
                    elif name in model_values and t in model_values[name]:
                        val = model_values[name][t]
                        ss_val = steady_states.get(name, 0.0)
                        e[j] = val - ss_val if not np.isnan(ss_val) else val

            curr = X @ curr + Y @ e
            trajectories[draw_idx, t, :] = curr

    x0 = dr.x0 if dr.x0 is not None else np.zeros(n_endo, dtype=float)
    return RandomSimulation(
        data=trajectories,
        variables=endo_names,
        steady_state=x0,
        units=units,
        canonical_units="deviation",
    )

def sim_to_nsim(sim: Any, units: UnitsType | None = None) -> pd.DataFrame:
    """Convert any simulation object, dict of DataFrames, or DataFrame to tidy long format
    with columns ``['shock', 't', 'variable', 'value']``.
    """
    import pandas as pd

    if isinstance(sim, SimulationResult):
        sim_dict = sim.to_dict(units=units)
    elif isinstance(sim, dict):
        sim_dict = sim
    elif isinstance(sim, pd.DataFrame):
        sim_dict = {"simulation": sim}
    else:
        raise TypeError(f"Unsupported simulation type for sim_to_nsim: {type(sim)}")

    cleaned: dict[Any, pd.DataFrame] = {}
    for k, df in sim_dict.items():
        if isinstance(df, pd.Series):
            df = df.to_frame()
        if "t" in df.columns:
            cleaned[str(k)] = df.drop(columns=["t"])
        else:
            cleaned[str(k)] = df

    if not cleaned:
        return pd.DataFrame(columns=["shock", "t", "variable", "value"])

    pdf = pd.concat(cleaned).reset_index()
    ppdf = pdf.rename(columns={"level_0": "shock", "level_1": "t"})
    ppdf = ppdf.melt(id_vars=["shock", "t"])
    return ppdf