dsgelab

Why it exists

The package is built so that the stochastic block and the deterministic transition are the same model.

In applied work, a perturbation solution is often built on one set of equations while the transition path is solved by a hand-written fixed point on another set, similar to the first, written by the same person, and not identical. Here both call the same residual function, so the equations behind them cannot differ.

Install

pip install dsgelab
pip install "dsgelab[symbolic]"

The package depends only on numpy and scipy. sympy is optional and serves the text front end, which supplies exact analytic derivatives.

dsgelab --example rbc        # solves a bundled model and prints the diagnostics

Examples

The examples run as written, in the order they appear, right after pip install. A checker (tools/check_site_examples.py) extracts every block on this page and runs it, so a block that stops running makes the check fail. Papers rerun in the package, with the number the page checks, are in Replications.

1. Solve a model that ships with the package

import dsgelab as dsg
from dsgelab.examples import rbc

model, params, guess = rbc.build()
sol = dsg.solve(model, params, guess=guess)

print(sol.summary())
irf = sol.irf("eps_a", size=0.01, periods=40)     # per cent deviations
mom = sol.moments(cov={"eps_a": 0.007})           # std, correlations, autocorrelations
dec = sol.variance_decomposition(cov={"eps_a": 0.007})

2. Declare your own model

A model consists of three lists of names and a function. The function reads variables by name, from cur for period t and from nxt for t+1.

import numpy as np
from dataclasses import dataclass
from dsgelab import Model, solve

@dataclass
class P:
    alpha: float = 0.36
    beta:  float = 0.99
    delta: float = 1.0
    gamma: float = 1.0
    rho:   float = 0.95

def equations(nxt, cur, eps, p):
    y  = cur.a * cur.k ** p.alpha
    y1 = nxt.a * nxt.k ** p.alpha
    return [
        nxt.k - ((1 - p.delta) * cur.k + y - cur.c),                    # resource
        cur.c ** -p.gamma - p.beta * nxt.c ** -p.gamma                  # Euler
            * (p.alpha * y1 / nxt.k + 1 - p.delta),
        np.log(nxt.a) - p.rho * np.log(cur.a) - eps.eps_a,              # technology
    ]

mine = Model(states=["k", "a"], jumps=["c"], shocks=["eps_a"], equations=equations)
sol_mine = solve(mine, P(), guess={"k": 0.2, "a": 1.0, "c": 0.36})

sol_mine.P                    # transition of the states
sol_mine.accuracy()           # the implied order must come out ~2

When variables are read from a vector by position, swapping one index for another raises no error, and the program silently solves a different model from the one intended. Reading each variable by name rules that mistake out.

3. The same equations, as text, with exact derivatives

from dsgelab.symbolic import from_equations

symbolic = from_equations(
    [
        "k(+1) = (1-delta)*k + a*k**alpha - c",
        "c**(-gamma) = beta * c(+1)**(-gamma) * (alpha*a(+1)*k(+1)**(alpha-1) + 1 - delta)",
        "log(a(+1)) = rho*log(a) + eps_a",
    ],
    states=["k", "a"], jumps=["c"], shocks=["eps_a"],
    parameters=["alpha", "beta", "rho", "delta", "gamma"],
)
sol_sym = dsg.solve(symbolic, P(), guess={"k": 0.2, "a": 1.0, "c": 0.36})

dsg.compare_irfs(sol_mine, sol_sym, "eps_a")["worst"]   # the two are the same model

The residual, the Jacobian and the Hessian are differentiated exactly. Exact derivatives matter at second order because a second difference divides by the square of the step, so a numerical Hessian carries about half the significant digits of the residual it came from.

4. Second order and risk corrections

so = dsg.solve_second_order(sol, cov={"eps_a": 0.007})

so.g_ss, so.h_ss          # what first order sets to zero
so.mean_correction()      # ergodic mean, in closed form
so.ergodic_mean()         # the same thing, in levels
so.accuracy()             # the implied order must come out ~3
so.simulate(periods=5000, seed=1)   # pruned, by default

A first-order simulation keeps the mean of the states at the steady state by construction. When the shift in the decision rules is large, it can get the sign of a mean effect wrong as well as its magnitude, and welfare rankings depend on that sign.

5. Perfect-foresight transition, with announced policy

from dataclasses import replace

p_old, p_new = P(), replace(P(), beta=0.985)
ss0 = dsg.steady_state(mine, {"k": 0.2, "a": 1.0, "c": 0.36}, p_old)
ss1 = dsg.steady_state(mine, ss0.values, p_new)

T = 80
surprise = dsg.perfect_foresight(mine, ss0, ss1, params=p_new, T=T)

schedule = [p_old] * 10 + [p_new] * (T + 1 - 10)           # announced at t=0
announced = dsg.perfect_foresight(mine, ss0, ss1, params=schedule, T=T)

surprise["k"], announced["k"]              # levels, (T+1,)
announced["c"][0] != surprise["c"][0]      # the announcement shows up on impact

The transition is solved by a stacked-time Newton method on the same residual function that the perturbation linearises. A surprise is a change in terminal with constant parameters, and a statutory schedule known in advance is passed as a sequence of parameter values, one per period.

6. Calibrate parameters on targets, and know whether they are identified

def observables(ss, p):
    y = ss["a"] * ss["k"] ** p.alpha
    return {"ky": ss["k"] / y, "iy": p.delta * ss["k"] / y}   # capital-output, investment share

res = dsg.calibrate(
    mine, P(),
    free={"beta": 0.95, "delta": 0.5},
    targets={"ky": 10.0, "iy": 0.25},
    guess={"k": 0.2, "a": 1.0, "c": 0.36},
    observables=observables,
    bounds={"beta": (0.5, 0.999), "delta": (0.001, 1.0)},
)
print(res.report())

ident = dsg.target_jacobian(mine, replace(P(), **res.values), free=["beta", "delta"],
                            targets=["ky", "iy"], guess={"k": 0.2, "a": 1.0, "c": 0.36},
                            observables=observables)
ident["identified"]        # True: two targets pin the two parameters

one = dsg.target_jacobian(mine, P(), free=["beta", "alpha"], targets=["k"],
                          guess={"k": 0.2, "a": 1.0, "c": 0.36})
one["identified"]          # False: two parameters, one target
one["null_directions"]     # the direction no target sees

The capital-output ratio and the investment share set β and δ together. In closed form δ = iy/ky = 0.025 and β = 1/(α/ky + 1 − δ) = 0.98912, which is what calibrate returns. Evaluated at that solution, target_jacobian reports that the two targets identify the two parameters. With β and α free and the capital stock as the only target, the elasticity matrix has rank one, and null_directions gives the combination of the two parameters that leaves the target unchanged. Done by hand, a calibration becomes nested bisections, one per parameter, each treating its target as if one parameter set it, although the targets are jointly determined. When that matrix lacks full rank, the calibration is not identified, and tightening the solver tolerance does not change that.

7. Independent checks

dsg.blanchard_kahn_report(sol)["margin"]        # the margin, not just "passed"

dsg.check_slack(sol, {"c_positive": lambda v: v["c"]}, size=-0.5)   # along the IRFs

dsg.steady_state_invariance(                    # an invariance that is a theorem of the model
    lambda v: dsg.steady_state(mine, {"k": 0.2, "a": 1.0, "c": 0.36}, replace(P(), rho=v)),
    [0.5, 0.8, 0.95], ["k", "c"])["ok"]         # k does not depend on the persistence

8. Derive the first-order conditions (FOCs) instead of writing them

You declare only the problem, and the package differentiates the Lagrangian. Because the derivatives are taken symbolically, the sign errors and omitted terms of differentiation by hand cannot enter the conditions. Those errors do not show up in the steady state, so a steady state that solves does not rule them out.

from dsgelab.blocks import Problem, optimality_conditions

prob = Problem(
    "log(c)",                                   # period return
    controls=["c", "i"], states=["k"], exogenous=["a"],
    constraints=["c + i = a*k**alpha",          # resources
                 "k(+1) = (1-delta)*k + i"],    # accumulation
    discount="beta", parameters=["alpha", "beta", "delta", "rho"],
)
foc = optimality_conditions(prob)

for eq in foc["equations"]:
    print(f"{eq:60s}  # {foc['detail'][eq]}")   # each equation says where it came from

derived = dsg.from_equations(
    foc["equations"] + ["log(a(+1)) = rho*log(a) + eps_a"],
    states=["k", "a"], jumps=foc["jumps"], shocks=["eps_a"],
    parameters=["alpha", "beta", "delta", "rho"])

Transversality conditions, inequality constraints and problems whose Euler equation does not follow from the derivative of next period's Lagrangian with respect to the inherited state are left to the modeller. The package derives the conditions of the problem that was declared, so an error in the declaration passes into the equations. Choosing the problem remains the modeller's task.

9. Third order

to = dsg.solve_third_order(so, cov={"eps_a": 0.007})

to.g_kkk, to.h_kkk        # the cubic terms
to.accuracy()             # the implied order must come out ~4
to.simulate(periods=5000, seed=1)    # pruned in three layers

10. Export to Dynare and to MacroModelling.jl

ss_sym = dsg.steady_state(symbolic, {"k": 0.2, "a": 1.0, "c": 0.36}, P())

text = dsg.to_dynare(symbolic, P(), steady=ss_sym, cov={"eps_a": 0.007},
                    path="model.mod", order=2)
jl = dsg.to_macromodelling(symbolic, P(), cov={"eps_a": 0.007}, path="model.jl")

"predetermined_variables k;" in text      # the timing is translated automatically

In dsgelab the state at t is the one available at t. In Dynare the predetermined variable is the one decided at t. A conversion by hand has to translate one timing into the other. The exporter does this translation, and it also writes exogenous processes one period back, so that the innovation appears without a lag. In the round trip the package writes the .mod and Dynare solves it. On the three reference models the two solutions differ by less than 1e-13 in steady state, IRFs and moments.

11. Optimal policy: the Ramsey planner

nk = dsg.ramsey_model(
    ["x - x(+1) + (1/sigma)*(i - pi(+1)) = 0",      # IS
     "pi = beta*pi(+1) + kappa*x + u",              # Phillips
     "u(+1) = rho*u + e"],                          # cost shock: exogenous process
    states=["u"], jumps=["x", "pi", "i"], shocks=["e"],
    parameters=["beta", "sigma", "kappa", "lam", "rho"],
    instruments=["i"],                              # the variable with no equation of its own
    objective="-(pi**2 + lam*x**2)/2", discount="beta")
pn = dict(beta=0.99, sigma=1.0, kappa=0.1, lam=0.25, rho=0.5)
ss_nk = dsg.ramsey_steady_state(nk, dict(u=0, x=0, pi=0, i=0), pn)
r = dsg.solve(nk, pn, steady=ss_nk).irf("e", size=1.0, percent=False)

x_lag = np.r_[0.0, r["x"][:-1]]
abs(r["pi"] + pn["lam"] / pn["kappa"] * (r["x"] - x_lag)).max() < 1e-12   # targeting rule

The planner's conditions are derived symbolically from the private-sector equations, and the result is an ordinary Model, with the multipliers and their lags as variables. The perspective is the timeless one, as in Dynare's ramsey_model. In this New Keynesian example the response follows Galí's targeting rule, π_t = −(λ/κ)(x_t − x_{t−1}), with an error below 1e-12.

12. Discretion, optimal simple rule and welfare

nk_eq = ["x - x(+1) + (1/sigma)*(i - pi(+1)) = 0",
         "pi = beta*pi(+1) + kappa*x + u",
         "u(+1) = rho*u + e"]
nk_kw = dict(states=["u"], jumps=["x", "pi", "i"], shocks=["e"],
             parameters=["beta", "sigma", "kappa", "lam", "rho"], instruments=["i"],
             objective="-(pi**2 + lam*x**2)/2", discount="beta")

disc = dsg.discretionary_policy(nk_eq, params=pn, cov={"e": 0.01}, **nk_kw)
w_disc = disc.welfare()["unconditional"]                       # without commitment
w_comm = dsg.welfare(nk, pn, steady=ss_nk, cov={"e": 0.01}).unconditional
w_disc < w_comm < 0                                            # commitment is worth something

rule = dsg.from_equations(["x - x(+1) + (1/sigma)*(i - pi(+1)) = 0",
                          "pi = beta*pi(+1) + kappa*x + u",
                          "i = phi_pi*pi + phi_u*u",            # the simple rule
                          "u(+1) = rho*u + e"],
                         states=["u"], jumps=["x", "pi", "i"], shocks=["e"],
                         parameters=["beta", "sigma", "kappa", "rho", "phi_pi", "phi_u"])
ot = dsg.osr(rule, dict(pn, phi_pi=1.5, phi_u=0.0), rule_params=["phi_pi", "phi_u"],
            weights={"pi": 1.0, "x": 0.25}, cov={"e": 0.01},
            guess=dict(u=0, x=0, pi=0, i=0))
ot.loss < ot.initial_loss

Discretion is solved by the algorithm of Dennis (2007), which Dynare's discretionary_policy also uses, and on the same model it reproduces Galí's closed-form solution. As in osr, the simple rule minimises the quadratic loss computed from the unconditional moments. Second-order welfare, conditional and unconditional, is defined as in evaluate_planner_objective.

13. News and surprises under perfect foresight

rbc_pf = dsg.from_equations(
    ["1/c = beta/c(+1)*(alpha*a(+1)*k(+1)**(alpha-1) + 1 - delta)",
     "c + k(+1) = a*k**alpha + (1-delta)*k",
     "a(+1) = 1 - rho + rho*a + e"],
    states=["k", "a"], jumps=["c"], shocks=["e"],
    parameters=["alpha", "beta", "delta", "rho"])
ppf = dict(alpha=0.33, beta=0.99, delta=0.025, rho=0.9)
kss = ((1/0.99 - 1 + 0.025)/0.33)**(1/(0.33 - 1))
ss_pf = dsg.steady_state(rbc_pf, dict(k=kss, a=1.0, c=kss**0.33 - 0.025*kss), ppf)

shock_path = np.zeros((80, 1)); shock_path[9, 0] = 0.02       # known from the start
anticipated = dsg.perfect_foresight(rbc_pf, ss_pf, ss_pf, params=ppf, T=80, shocks=shock_path)
surprise = dsg.perfect_foresight_with_expectation_errors(       # learned only at t = 10
    rbc_pf, ss_pf, ss_pf, params=ppf, T=80,
    news=[dsg.News(period=10, shocks=shock_path)])

anticipated["c"][0] != ss_pf.values["c"]     # consumption reacts to the announcement
abs(surprise["c"][0] - ss_pf.values["c"]) < 1e-12   # and not to what nobody knows

A news item is a new expected path of shocks and parameters, and a new terminal steady state. At each news item the agents replan from the states already realised. The procedure is that of Dynare's perfect_foresight_with_expectation_errors_solver, and a .mod with the information table in CSV becomes the list of news items through ModFile.news.

14. Bayesian estimation

ar = dsg.from_equations(["a(+1) = rho*a + e", "y = 1 + a"], states=["a"], jumps=["y"],
                       shocks=["e"], parameters=["rho"])
sol_ar = dsg.solve(ar, dict(rho=0.8), guess=dict(a=0, y=1))
series = sol_ar.to_levels(sol_ar.simulate(periods=400, cov={"e": 0.1}, seed=1))["y"]

pb = dsg.EstimationProblem(
    model=ar, data={"y": series}, params=dict(rho=0.5), observables=["y"],
    priors={"rho": ("beta", 0.5, 0.2), "stderr e": ("inv_gamma", 0.1, float("inf"))},
    guess=dict(a=0, y=1))
mode = dsg.posterior_mode(pb)                       # mode, Hessian, Laplace
chain = dsg.metropolis_hastings(pb, mode, draws=1500, chains=2, seed=1)

abs(mode.values["rho"] - 0.8) < 0.1
max(chain.rhat().values()) < 1.05
abs(chain.log_marginal_mhm() - mode.laplace) < 0.5

EstimationProblem holds the model, the data, the priors and the observables. The priors are Dynare's, set by mean and standard deviation. Observables may be in levels or in logs, may carry measurement error and may have missing values. From the problem the package computes the mode with the Laplace approximation, maximum likelihood, Metropolis with R-hat and effective sample size, SMC, the prior predictive check and the marginal likelihood. ModFile.estimation_problem builds an EstimationProblem from a .mod with estimated_params and varobs.

15. Occasionally binding constraints (OccBin)

credit = dsg.OccbinModel.from_equations(
    ["c = y + b - R*LB",
     "LB(+1) = b",
     "[name='debt', relax='slack'] b = m*y",        # the limit holds ...
     "[name='debt', bind='slack'] lam = 0",         # ... or the multiplier is zero
     "lam = 1/c - beta*R/c(+1)",
     "log(y(+1)) = rho*log(y) + u"],
    constraints={"slack": {"bind": "lam < 0", "relax": "b > m*y"}},
    states=["LB", "y"], jumps=["c", "b", "lam"], shocks=["u"],
    parameters=["R", "beta", "m", "rho"])
pc = dict(R=1.04, beta=0.95, m=0.8, rho=0.85)
c0 = 1 + 0.8 - 1.04*0.8
ob = dsg.occbin_simulate(credit, pc, {"u": {2: 0.06}}, periods=40,
                        guess=dict(LB=0.8, y=1.0, b=0.8, c=c0, lam=(1 - 0.95*1.04)/c0))

ob.binding["slack"].any()                                   # high income loosens the limit
max(abs(ob.piecewise["lam"][ob.binding["slack"]])) < 1e-14  # and the multiplier goes to zero

Each regime is linearised at the reference steady state, with its residual as a constant. For a guess of the dates at which each constraint holds, the path follows time-varying rules found backwards, and the guess is then revised by the rule of Dynare's +occbin, with curb_retrench if requested. A .mod with tagged equations (bind=, relax=), occbin_constraints and shocks(surprise) runs through ModFile.occbin_simulate().

16. Extended path, deterministic and stochastic

tree = dsg.from_equations(                         # Burnside (1998): the price of a Lucas tree
    ["x(+1) = (1 - rho)*xbar + rho*x + e",
     "1 = beta*exp(theta*x(+1))*(1 + y(+1))/y"],
    states=["x"], jumps=["y"], shocks=["e"], parameters=["beta", "theta", "rho", "xbar"])
pa = dict(beta=0.95, theta=-1.5, rho=-0.139, xbar=0.0179)
g0 = 0.95*np.exp(-1.5*0.0179)
ss_a = dsg.steady_state(tree, dict(x=0.0179, y=g0/(1 - g0)), pa)
draws = 0.0348*np.random.default_rng(0).standard_normal((20, 1))

ep = dsg.extended_path(tree, ss_a, pa, draws, horizon=60)                     # no risk ahead
sep = dsg.extended_path(tree, ss_a, pa, draws, horizon=60, order=1,           # next-period risk
                       cov=[[0.0348**2]])
bool((sep["y"][1:] > ep["y"][1:]).all())       # with theta < 0, risk raises the tree's price

Each period the shock arrives, the states it leads to are computed, and the agents solve a nonlinear problem that reaches the steady state horizon periods later. With order > 0 the shocks of the following periods are integrated by Gauss–Hermite quadrature over a tree of plans, sparse or full, with the optional second-order risk correction (hybrid). The algorithm and these options are those of Dynare's extended_path.

17. Declaring the model: lags, calibration inside the model, homotopy, diagnostics

rbc_cal = dsg.from_equations(
    ["1/c = beta/c(+1)*(alpha*k(+1)**(alpha-1) + 1 - delta)",
     "c + k(+1) = k**alpha + (1-delta)*k",
     "q = k**alpha",
     "dq = q - q(-1)"],                         # the lag of a jump becomes an auxiliary state
    states=["k"], jumps=["c", "q", "dq"], shocks=[], parameters=["alpha", "beta", "delta"],
    calibration_equations={"beta": "k/q = 10"})  # beta is what sets k/q to 10
ss_cal = dsg.steady_state(rbc_cal, dict(k=10, c=1, q=1.5, dq=0),
                         dict(alpha=0.33, beta=0.99, delta=0.025))
rbc_cal.auxiliary                               # {'LAG1_q': ('lag', 'q')}
abs(ss_cal.calibrated["beta"] - 1/(0.33/10 + 1 - 0.025)) < 1e-12
dsg.diagnose(rbc_cal, ss_cal.params, ss_cal.values)["ok"]

rbc_h = dsg.from_equations(
    ["1/c = beta/c(+1)*(alpha*k(+1)**(alpha-1) + 1 - delta)",
     "c + k(+1) = k**alpha + (1-delta)*k"],
    states=["k"], jumps=["c"], shocks=[], parameters=["alpha", "beta", "delta"])
ss_h = dsg.steady_state(rbc_h, dict(k=30, c=2), dict(alpha=0.33, beta=0.99, delta=0.9),
                       homotopy={"delta": (0.02, 0.9)})   # the guess is good for delta = 0.02

A closed-form steady state (steady_state_model, here or in a .mod) is checked against the equations before it is used and, as in Dynare, refused if it does not set them to zero. dsg.assemble combines the problems of several agents and the market-clearing conditions into a single model.

18. Trends and balanced growth

trended = dsg.from_mod("""
trend_var(growth_factor = g) A;
var(deflator = A) c k y;
parameters alpha beta delta g;
alpha = 0.33; beta = 0.99; delta = 0.025; g = 1.005;
model;
1/c = beta/c(+1)*(alpha*y(+1)/k + 1 - delta);
c + k = y + (1-delta)*k(-1);
y = k(-1)^alpha*A^(1-alpha);
end;
initval; c = 0.8; k = 10; y = 1; end;
""")
ss_t = trended.steady_state()                    # the model comes out already detrended

solow = dsg.from_equations(
    ["E(+1) = gE*E", "y = k**alpha*E**(1-alpha)", "k(+1) = (1-delta)*k + s*y"],
    states=["E", "k"], jumps=["y"], shocks=[], parameters=["alpha", "delta", "s", "gE"])
bgp = dsg.balanced_growth(solow, dict(alpha=0.33, delta=0.02, s=0.2, gE=1.01))
abs(bgp.growth["y"] - 1.01) < 1e-12
abs(bgp.ratio("k", "y") - 0.2/(1.01 - 1 + 0.02)) < 1e-10

In a .mod written in levels, trend_var and var(deflator=) are detrended as Dynare does it, and an equation that still depends on the level of the trend is refused. balanced_growth computes the growth factors and the ratios between variables that share the trend, which is the problem Dynare's bgp.write solves.

What is in the box

The one convention you have to get right

eps is the innovation dated t+1, unknown at t, with E_t eps = 0. An AR(1) is written as

nxt.a - p.rho * cur.a - eps.eps_a

Validation

The validation has three layers. The test suite compares the package with closed forms and answers known analytically. The scripts in validation/ compare it with Dynare on Dynare's own examples and test files, with MacroModelling.jl on one calibration equation, and with the original MATLAB code of Heutel (2012), shown in Replications. A last comparison, on an applied model from an unpublished paper by the author, cannot be rerun from the repository.

In the second layer both the models and the reference software come from other authors. The third uses the same reference software on a model of the author's own. Each feature is validated against an independent reference where one exists (another implementation, a closed form or a published result). A few rest on internal tests only, namely autocorrelations, variance decomposition, pruning, the second-order generalised IRF, announced parameter paths and the Blanchard–Kahn check, and the export to MacroModelling.jl has no round trip.

Against closed forms

Of the package's 148 tests, those in this layer compare its output with answers that can be derived on paper. With log utility and full depreciation, the RBC solution must equal P = [[α, 1], [0, ρ]]. In the three-equation New Keynesian model with φx = 0.125, the Taylor rule must give a determinate equilibrium at φπ = 1.5 and an indeterminate one at φπ = 0.8. On a lognormal AR(1), the second-order coefficient h_ss has to equal sigma², so the decision rule shifts by sigma²/2, the Jensen term. The perfect-foresight transition is compared with the closed-form path k' = αβk^α. The order-of-accuracy test requires that halving the displacement divide the first-order residual by four and the second-order residual by eight.

Against Dynare's reference models

The examples from dynare.org/resources were run in Dynare and compared with dsgelab variable by variable.

Model	Steady state	IRFs	Theoretical moments
example1 (Collard, correlated shocks)	1.2e-14	2.9e-14	2.6e-14	OK
example2 (the same in logarithms)	2.1e-05	4.3e-06	3.3e-06	OK*
RBC_baseline (Pfeifer, 15 variables)	0.0e+00	2.3e-14	3.5e-14	OK

* In example2 Dynare's steady stops early, leaving a residual of 1.05e-06 in the model equations against 1.03e-44 for dsgelab, and the whole discrepancy follows from that residual. Because example1 and example2 are the same economy in different coordinates, they allow an independent check, and exp(s.s. of example2) matches example1 in levels to 1e-15.

Dynare's examples, read without porting

dsg.from_mod reads the .mod as it is (macros, local variables, leads and lags of any order, shocks inside any equation). The steady state, IRFs and theoretical standard deviations are then compared with what Dynare obtains from the same file.

Model	Variables	Steady state	IRFs	Moments
fs2000 (Schorfheide 2000)	14	9.2e-15	3.4e-13	3.4e-13	OK
NK_baseline (Fernández-Villaverde and Rubio-Ramírez)	28	9.3e-14	5.9e-11	2.2e-11	OK
Gali_2015 (ch. 3, unit root)	14	0.0e+00	1.6e-14	7.5e-15	OK
example3 (Collard)	6	8.3e-15	1.4e-14	2.2e-14	OK
agtrend (Aguiar–Gopinath, macros)	15	3.9e-08	2.1e-07	2.3e-07	OK*
fs2000_nonstationary (in levels, trend_var and deflator)	11	1.4e-08	1.5e-08	2.2e-08	OK*

* The agreement is limited by the residual of 4e-07 that Dynare's steady state leaves in the equations (7e-08 in fs2000_nonstationary, which Dynare and dsgelab each detrend on their own). The three models of the previous table also go through this route. For example1 and RBC_baseline the deviations stay below 1.2e-13, the largest being the steady state of RBC_baseline, 1.1e-13, and example2 gives the numbers of that table.

Optimal policy against Dynare's ramsey_model

The comparison covers the steady state and IRFs of the private variables and of the multipliers, which both programs number by the same equation and with the same sign.

Model	Steady state	IRFs
nk_linear (Galí, ch. 5)	0.0e+00	5.3e-15	OK
rbc_labortax (capital as a chosen stock)	5.0e-11	4.9e-13	OK
nk_ramsey_model (Dynare's tests, nonlinear)	2.7e-12	2.1e-14	OK

Discretion, simple rules and welfare against Dynare

from_mod reads each file without changes. Three come from Dynare's tests, and nk_osr_planner.mod and rbc_labortax.mod were written in this repository.

Test	What is compared	Deviation
osr_example.mod	loss at Dynare's optimum	2.6e-16	OK*
nk_osr_planner.mod	optimal simple rule on welfare	1.5e-16	OK
dennis_1.mod	IRFs under discretion; unconditional and conditional welfare	4.3e-08; 5.3e-10, 6.1e-10	OK
nk_ramsey_model.mod	Ramsey read from the .mod; 2nd-order welfare	2.6e-14; 1.4e-12	OK
rbc_labortax.mod	2nd-order welfare with capital	3.5e-13	OK

* The dsgelab optimum has a lower loss than Dynare's, 0.4266626 against 0.4268336. Evaluated at the dsgelab point, Dynare's loss differs from dsgelab's by 4.3e-15, so the difference comes from Dynare's optimiser stopping early on a flat surface. At the optimum the rule neutralises the demand shock, with gamma_y_ = 1/|σ| = 16.667. The 4e-08 floor in the discretion IRFs is Dynare's tolerance, sqrt(eps).

Perfect foresight against Dynare

Each file in Dynare's tests/deterministic_simulations is read by from_mod, and the comparison covers every variable in every period.

Test	What it exercises	Deviation
rbc_det1.mod	initial condition (histval)	1.1e-14	OK
rbc_det2.mod	shock in the first period	6.6e-15	OK
rbc_det3.mod	anticipated shocks, periods 4 to 8	8.0e-14	OK
pfwee.mod	expectation errors, transitory and permanent shocks	2.8e-12	OK
pfwee_constant_sim_length.mod	the same, constant horizon	3.0e-12	OK

By default Dynare stops at a residual of 1e-5 and recomputes the terminal steady state with tolerance eps^(1/3). With those defaults, the terminal k of pfwee misses the closed-form solution by 4e-5. The validation runs Dynare at 1e-12.

Estimation against Dynare

The model is fs2000.mod (Schorfheide, 2000) from Dynare's examples/ folder, read by from_mod with its priors and data.

What is compared	Deviation
likelihood and prior density, at the initial point	6e-12; 4e-13	OK
smoothed variables; smoothed shocks	8.4e-15; 1.1e-15	OK
shock_decomposition and realtime_shock_decomposition	5e-15	OK
8-period forecast; 90% bands	2e-15; 6e-09	OK
posterior mode (dsgelab kernel ≥ Dynare's)	+6e-06	OK
Laplace with Dynare's covariance	2e-13	OK
Metropolis means, in units of Monte Carlo error (largest 2.89)	< 3	OK
diffuse filter (observed price level), at two points	2.5e-12	OK
conditional forecast (tests/conditional_forecasts)	1.6e-13	OK

The likelihood matches to the last digit only with Dynare's riccati_tol switched off. By default Dynare switches to the filter's steady-state gain when the variance converges to 1e-6, which moves the number in the seventh decimal place. The Laplace approximation with dsgelab's own Hessian differs from Dynare's by 0.03, and that difference comes from the two numerical Hessians.

Nonlinearities and model declaration against Dynare

The files in tests/occbin, tests/ep, tests/bgp and tests/trend_var of Dynare 7.2 are read by from_mod, and each variable is checked in each period.

Test	What is compared	Deviation
borrcon.mod	OccBin, one constraint, curb_retrench: piecewise and linear paths; regimes	2.2e-16; 30/30	OK
dynrbc.mod	OccBin, two constraints and four regimes, 52 surprise shocks	1.6e-15; 169/169	OK
burnside.mod	extended path; stochastic of order 1 and 2, sparse and full trees, hybrid	1.1e-11	OK*
linearmodel0.mod	extended path on a linear model; Dynare's first-order solution	2.0e-15; 3.1e-14	OK
solow, ramsey, nk	the bgp.write problem: residuals at five points; the tests' assertions	2.2e-16; 2.8e-12	OK
walras.mod, bad_steady.mod, homotopy.mod	model_diagnostics; steady_state_model refused; homotopy	same report	OK

* Dynare is run here with a solver tolerance of 1e-12, against the default of 1e-5, and that tolerance is what limits the agreement.

A calibration equation against MacroModelling.jl

In a small RBC model the steady-state target k/q = 2.5 sets β, and the model was solved by dsgelab and by MacroModelling.jl. Both return β = 0.847457627118644, the closed-form value, and their steady-state values differ by at most 4.4e-15.

Against Dynare and MacroModelling.jl, on an applied model by the author

In a DSGE model with 71 variables, three income groups and five shocks, the second-order mean corrections of dsgelab match Dynare to 1.1e-05 percentage points and MacroModelling.jl to 1.2e-05, at the October 2026 calibration. Dynare uses unpruned second order and MacroModelling uses pruned second order. dsgelab solves the Sylvester equation and returns the mean in closed form, without simulating. The model belongs to an unpublished paper and is not in this repository, so this comparison cannot be rerun from here, unlike the ones above.

How it compares with what already exists

Several Python packages solve DSGE models, and the table sets their capabilities side by side with those of dsgelab.

Capability	dsgelab	gEconpy	dolo	PyMacLab	econpizza
Derives the FOCs from the agent's problem	yes	yes	—	—	—
1st-order perturbation	yes	yes	yes	yes	—
2nd order + risk corrections	yes	no	yes	yes	—
3rd order (deterministic block)	yes	—	yes	—	—
3rd order with risk in the slope	yes (g_ssk, h_ssk); Dynare's full tensor, no	—	yes	—	—
Global / projection methods	partial (Chebyshev collocation, up to three states)	—	yes	—	partial
Heterogeneous agents	—	—	—	—	yes
Nonlinear perfect foresight	yes	yes	yes	?	yes
Bayesian estimation	yes (mode, Metropolis–Hastings and SMC, with the Kalman filter on the first-order solution)	yes	—	—	—
Calibration on targets + identification	yes	partial	?	?	?
Export to Dynare and MacroModelling.jl	yes	Dynare only	Dynare only	—	—
Self-verified order of accuracy	yes (orders 2, 3, 4)	?	?	?	?
Dependencies	numpy, scipy	sympy, PyMC	pandas (old pin)	numpy, scipy	JAX
Installs on Python 3.14	yes	yes	no	no	?

Survey of 9 September 2026, revised on 11 and 12 October 2026, with gEconpy 2.2. Later releases of these packages may change some cells, and docs/comparison.md gives the source of each one.

In dsgelab a declared problem goes from the derived first-order conditions to a solution and to an export for Dynare or MacroModelling.jl, with the timing translated automatically. The same residual function serves the perturbation and the deterministic transition. The second-order solution runs within Python, without calling Dynare through Octave. Blanchard–Kahn and the convergence of the steady state are checked on every solve, the order of accuracy and the identification of a calibration take one call each, and the only dependencies are numpy and scipy. In the three uses named in the note above (Bayesian estimation through PyMC, global methods and heterogeneous agents), the packages named there offer more than dsgelab. The full table, which separates the cells verified by installing a package from those read only in its documentation, is in docs/comparison.md.

Citing

If the package contributed to a published result, it belongs in the references, as any other methodological input does. The reference should include the version, because a number that reproduces only under one version cannot be reproduced without knowing which.

BibTeX

@software{dsgelab_2026,
  author  = {Moreira, Gustavo Hoffmann},
  title   = {{dsgelab}: a general toolkit for solving DSGE models in Python},
  year    = {2026},
  version = {0.5.0},
  license = {MIT},
  doi     = {10.5281/zenodo.23456475},
  url     = {https://github.com/GustavoHM/dsgelab},
  note    = {Python package}
}

APA

Moreira, G. H. (2026). dsgelab: A general toolkit for solving DSGE models in Python
(Version 0.5.0) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.23456475

In the text

Next to the result, give the version along with the name of the package, as in the sentence below.

The model was solved with dsgelab 0.5.0 (Moreira, 2026), by second-order perturbation.

Citing the tools this package uses

Cite Dynare when the model was exported to it or checked there, and MacroModelling.jl in the same case. NumPy and SciPy ask explicitly to be cited, and every computation in dsgelab runs on them.

Documentation