DSGE Models

MacroEconometricModels.jl provides a complete toolkit for specifying, solving, simulating, and estimating Dynamic Stochastic General Equilibrium (DSGE) models. The package covers the full workflow from model definition through structural estimation, with seven solution methods spanning linear, higher-order, and global approaches.

• Specification: The @dsge macro provides a domain-specific language for writing equilibrium conditions with time-indexed variables
• Steady State: Analytical or numerical steady-state computation via NonlinearSolve.jl (TrustRegion() default) with built-in constrained solvers (Optim.jl, NLopt.jl) and optional JuMP backends
• Linearization: Automatic first-order approximation via numerical Jacobians in the Sims (2002) canonical form
• Linear Solvers: Three first-order solvers –- Gensys (Sims 2002), Blanchard-Kahn (1980), and Klein (2000) –- producing the state-space solution; see Linear Solvers
• Nonlinear Methods: Up to 3rd-order perturbation with Andreasen, Fernandez-Villaverde & Rubio-Ramirez (2018) pruning, Chebyshev collocation, policy function iteration, and value function iteration (with Howard steps and Anderson acceleration) for globally accurate policy functions; see Nonlinear Methods
• Constraints: Perfect foresight paths, OccBin occasionally-binding constraints (Guerrieri & Iacoviello 2015), and constrained optimization via Optim.jl/NLopt.jl (built-in) with optional JuMP/Ipopt (NLP) and PATH (MCP) backends; see Constraints
• Estimation: Four GMM-based methods (one-step, two-step, iterative, CU) for IRF matching, plus Bayesian estimation via SMC, SMC$^2$ with two-stage delayed acceptance, and Random-Walk Metropolis-Hastings; see Estimation
• Simulation and IRFs: Stochastic and pruned simulation, analytical and generalized impulse responses, FEVD, and unconditional moments via Lyapunov equation; see Nonlinear Methods
• Historical Decomposition: Kalman smoother-based shock attribution for linear models, FFBSi particle smoother for nonlinear models, and Bayesian posterior bands; see Historical Decomposition

All results integrate with plot_result() for interactive D3.js visualization and report() for publication-quality output.

using MacroEconometricModels, Random
Random.seed!(42)
# Pre-define the RBC spec for reuse across blocks
_spec_rbc = @dsge begin
    parameters: β = 0.99, α = 0.36, δ = 0.025, ρ = 0.9, σ = 0.01
    endogenous: Y, C, K, A
    exogenous: ε_A
    Y[t] = A[t] * K[t-1]^α
    C[t] + K[t] = Y[t] + (1 - δ) * K[t-1]
    1 = β * (C[t] / C[t+1]) * (α * A[t+1] * K[t]^(α - 1) + 1 - δ)
    A[t] = ρ * A[t-1] + σ * ε_A[t]
    steady_state = begin
        A_ss = 1.0
        K_ss = (α * β / (1 - β * (1 - δ)))^(1 / (1 - α))
        Y_ss = K_ss^α
        C_ss = Y_ss - δ * K_ss
        [Y_ss, C_ss, K_ss, A_ss]
    end
end
_spec_rbc = compute_steady_state(_spec_rbc)

<< @setup-block not executed in draft mode >>

Quick Start

Recipe 1: Solve and plot IRFs

# Specify a simple RBC model
spec = @dsge begin
    parameters: β = 0.99, α = 0.36, δ = 0.025, ρ = 0.9, σ = 0.01
    endogenous: Y, C, K, A
    exogenous: ε_A

    Y[t] = A[t] * K[t-1]^α
    C[t] + K[t] = Y[t] + (1 - δ) * K[t-1]
    1 = β * (C[t] / C[t+1]) * (α * A[t+1] * K[t]^(α - 1) + 1 - δ)
    A[t] = ρ * A[t-1] + σ * ε_A[t]

    steady_state = begin
        A_ss = 1.0
        K_ss = (α * β / (1 - β * (1 - δ)))^(1 / (1 - α))
        Y_ss = K_ss^α
        C_ss = Y_ss - δ * K_ss
        [Y_ss, C_ss, K_ss, A_ss]
    end
end
spec = compute_steady_state(spec)

sol = solve(spec)
result = irf(sol, 40)

<< @example-block not executed in draft mode >>

plot_result(result)

Recipe 2: Second-order perturbation with pruning

psol = perturbation_solver(spec; order=2)
Y_sim = simulate(psol, 1000)  # pruned simulation (Kim et al. 2008)
girf = irf(psol, 40; irf_type=:girf, n_draws=100)
nothing # hide

<< @example-block not executed in draft mode >>

Recipe 3: Estimate with GMM

est = estimate_dsge(spec, Y_data, [:β, :α, :ρ];
                    method=:irf_matching, var_lags=4, irf_horizon=20)
report(est)

Recipe 4: OccBin with ZLB

constraint = parse_constraint(:(R[t] >= 0), spec)
occ_sol = occbin_solve(spec, constraint; shock_path=shocks)
occ_irf = occbin_irf(spec, constraint, 1, 40)
plot_result(occ_irf)

Recipe 5: Chebyshev projection

proj = collocation_solver(spec; degree=5, grid=:tensor, max_iter=200)
y = evaluate_policy(proj, proj.steady_state[proj.state_indices])
err = max_euler_error(proj)
nothing # hide

<< @example-block not executed in draft mode >>

Recipe 6: Bayesian estimation via SMC$^2$

using Distributions
result = estimate_dsge_bayes(spec, data, [0.99, 0.9, 0.01];
    priors=Dict(:β => Beta(99, 1), :ρ => Beta(5, 2), :σ => InverseGamma(2, 0.01)),
    method=:smc2, observables=[:y], n_smc=200, n_particles=100,
    solver=:projection, solver_kwargs=(degree=5,))
report(result)

Recipe 7: Historical decomposition

data_hd = simulate(sol, 100)
hd = historical_decomposition(sol, data_hd, [:Y, :C, :K, :A])
report(hd)

<< @example-block not executed in draft mode >>

Model Specification

The @dsge macro provides a domain-specific language for specifying DSGE models. It parses the model block into a DSGESpec{T} object containing equations, parameters, and variable declarations.

Syntax

spec_demo = @dsge begin
    parameters: β = 0.99, α = 0.36, δ = 0.025, ρ = 0.9, σ = 0.01
    endogenous: Y, C, K, A
    exogenous: ε_A

    # Equations (one per endogenous variable)
    Y[t] = A[t] * K[t-1]^α
    C[t] + K[t] = Y[t] + (1 - δ) * K[t-1]
    1 = β * (C[t] / C[t+1]) * (α * A[t+1] * K[t]^(α - 1) + 1 - δ)
    A[t] = ρ * A[t-1] + σ * ε_A[t]
end

<< @example-block not executed in draft mode >>

Blocks

Time Subscripts

Variables with [t+1] subscripts generate expectation errors in the Sims (2002) canonical form. The number of forward-looking equations determines the dimension of the $\Pi$ matrix and, via the Blanchard-Kahn (1980) condition, the number of unstable eigenvalues required for determinacy.

Return Value

Steady State

The steady state $\bar{y}$ satisfies the equilibrium system in the absence of shocks:

\[f(\bar{y}, \bar{y}, \bar{y}, 0, \theta) = 0\]

where:

• $\bar{y}$ is the $n \times 1$ vector of endogenous variables at the steady state
• $\theta$ is the vector of deep structural parameters
• $f$ is the system of $n$ equilibrium conditions

For the RBC model above, the analytical steady state is:

\[\bar{A} = 1, \quad \bar{K} = \left(\frac{\alpha\beta}{1 - \beta(1-\delta)}\right)^{\frac{1}{1-\alpha}}, \quad \bar{Y} = \bar{K}^\alpha, \quad \bar{C} = \bar{Y} - \delta\bar{K}\]

Numerical Computation

compute_steady_state uses NonlinearSolve.jl to solve the system $f(\bar{y}, \bar{y}, \bar{y}, 0, \theta) = 0$. The default algorithm is TrustRegion(), which is robust to poor starting points. Box constraints (e.g., non-negativity) are handled natively via NonlinearSolve's bounded problem formulation.

spec = compute_steady_state(spec)
report(spec)

<< @example-block not executed in draft mode >>

The solver converges to the steady state from a default initial guess of ones. For models with multiple equilibria, providing a good starting point via initial_guess avoids convergence to an economically irrelevant solution.

Analytical Steady State

For models where the steady state has a closed-form solution, specify it in a steady_state = begin ... end block. The block must return a vector matching the endogenous variable ordering:

spec = @dsge begin
    parameters: β = 0.99, α = 0.36, δ = 0.025, ρ = 0.9, σ = 0.01
    endogenous: Y, C, K, A
    exogenous: ε_A

    Y[t] = A[t] * K[t-1]^α
    C[t] + K[t] = Y[t] + (1 - δ) * K[t-1]
    1 = β * (C[t] / C[t+1]) * (α * A[t+1] * K[t]^(α - 1) + 1 - δ)
    A[t] = ρ * A[t-1] + σ * ε_A[t]

    steady_state = begin
        A_ss = 1.0
        K_ss = (α * β / (1 - β * (1 - δ)))^(1 / (1 - α))
        Y_ss = K_ss^α
        C_ss = Y_ss - δ * K_ss
        [Y_ss, C_ss, K_ss, A_ss]
    end
end
spec = compute_steady_state(spec)
nothing # hide

<< @example-block not executed in draft mode >>

When the steady_state block is provided, compute_steady_state (or solve) uses it directly and validates the result against the equations. The analytical path is faster and avoids numerical convergence issues, but the user is responsible for correctness –- the validator checks that $\|f(\bar{y})\| < 10^{-10}$.

Constrained Steady State

For models with variable bounds –- such as a zero lower bound on the nominal interest rate or non-negativity of consumption –- box constraints work out of the box. NonlinearSolve.jl attempts the solve first; if bounds are violated, the solver auto-escalates to Optim.jl Fminbox(LBFGS()):

# ZLB: interest rate cannot go below zero
bound = variable_bound(:R, lower=0.0)

# Solve constrained steady state (auto-escalates to Optim.jl if bounds bind)
spec = compute_steady_state(spec; constraints=[bound])

For nonlinear inequality constraints, NLopt.jl LD_SLSQP is the default solver –- no additional packages required:

# Debt-to-GDP ceiling: fn(y, ...) <= 0
debt_limit = nonlinear_constraint(
    (y, y_lag, y_lead, e, theta) -> y[debt_idx] / y[gdp_idx] - 0.6;
    label="Debt-to-GDP <= 60%"
)

spec = compute_steady_state(spec; constraints=[bound, debt_limit])

Advanced: Explicit Solver Backends

For large-scale problems or complementarity formulations, JuMP-based backends provide additional power:

# Ipopt (NLP): handles general nonlinear constraints
import JuMP, Ipopt
spec = compute_steady_state(spec; constraints=[bound], solver=:ipopt)

# PATH (MCP): natural for complementarity problems (e.g., ZLB)
import JuMP, PATHSolver
spec = compute_steady_state(spec; constraints=[bound], solver=:path)

Linearization

linearize computes a first-order Taylor expansion around the steady state using numerical Jacobians (central differences). It produces the Sims (2002) canonical form:

\[\Gamma_0 \, y_t = \Gamma_1 \, y_{t-1} + C + \Psi \, \varepsilon_t + \Pi \, \eta_t\]

where:

• $y_t$ is the $n \times 1$ vector of endogenous variables (deviations from steady state)
• $\Gamma_0$ is the $n \times n$ coefficient matrix on current-period variables
• $\Gamma_1$ is the $n \times n$ coefficient matrix on lagged variables
• $C$ is the $n \times 1$ constant vector
• $\Psi$ is the $n \times n_{shocks}$ shock loading matrix
• $\Pi$ is the $n \times n_{expect}$ expectation error selection matrix
• $\varepsilon_t$ is the vector of exogenous shocks
• $\eta_t = y_t - E_{t-1}[y_t]$ is the vector of expectation errors for forward-looking variables

ld = linearize(spec)
nothing # hide

<< @example-block not executed in draft mode >>

The matrix pair $(\Gamma_0, \Gamma_1)$ defines a generalized eigenvalue problem whose solution governs the model dynamics. The three Linear Solvers –- Gensys, Blanchard-Kahn, and Klein –- each decompose this pencil to extract the stable state-space representation.

Return Value

Complete Example

This example specifies, solves, and analyzes a full RBC model using the core functions covered on this page:

# Verify steady state (spec defined in Quick Start)
report(spec)

<< @example-block not executed in draft mode >>

# Linearize and inspect the canonical form matrices
ld = linearize(spec)

# Solve with default Gensys method
sol = solve(spec)

# IRFs and FEVD
result = irf(sol, 40)
report(result)

<< @example-block not executed in draft mode >>

plot_result(result)

The spec object stores the parsed model. linearize produces the Sims (2002) canonical form $(\Gamma_0, \Gamma_1, C, \Psi, \Pi)$. solve dispatches to the Gensys algorithm and returns a DSGESolution with the state-space representation $y_t = G_1 y_{t-1} + C + \text{impact} \cdot \varepsilon_t$. The IRF and FEVD functions operate on this solution. For higher-order or global solutions, see Nonlinear Methods.

Common Pitfalls

1. Steady-state validation failure: When providing an analytical steady_state block, the validator checks $\|f(\bar{y})\| < 10^{-10}$. A common cause of failure is mismatched variable ordering –- the returned vector must match the endogenous: declaration order exactly.

2. Equation count mismatch: The number of equations must equal the number of endogenous variables. Missing an equilibrium condition or double-counting a definition produces DimensionMismatch. Each equation is written as LHS = RHS; the parser rearranges to residual form automatically.

3. Timing convention confusion: $K_t$ chosen at time $t$ is written K[t]. Beginning-of-period capital (predetermined) is K[t-1]. A forward-looking Euler equation uses C[t+1]. Misplacing a time subscript silently changes the $\Gamma_0$/$\Gamma_1$ structure and can cause indeterminacy.

4. Numerical steady state converges to wrong equilibrium: For models with multiple equilibria, the default initial guess (vector of ones) may converge to an economically irrelevant solution. Provide initial_guess close to the desired equilibrium, or use the analytical steady_state block.

5. Constrained steady state: Box constraints (variable bounds) are handled by NonlinearSolve.jl without additional dependencies. Nonlinear inequality constraints require import JuMP, Ipopt. PATH MCP requires import JuMP, PATHSolver.

References

• Sims, C. A. (2002). Solving Linear Rational Expectations Models. Computational Economics, 20(1-2), 1-20. DOI

• Blanchard, O. J., & Kahn, C. M. (1980). The Solution of Linear Difference Models under Rational Expectations. Econometrica, 48(5), 1305-1311. DOI