dsgelab

Por que este pacote existe

O pacote foi projetado sob a restrição de que o bloco estocástico e a transição determinística sejam o mesmo modelo.

No trabalho aplicado, é comum que a solução por perturbação seja montada sobre um conjunto de equações e que a trajetória de transição seja resolvida, por um ponto fixo programado à mão, sobre outro conjunto, parecido com o primeiro e escrito pela mesma pessoa, mas diferente dele. Aqui as duas chamam a mesma função de resíduos e, por isso, não podem divergir.

Instalar

pip install dsgelab
pip install "dsgelab[symbolic]"

O pacote depende só de numpy e scipy. O sympy é opcional e serve ao front end textual, que calcula as derivadas analiticamente.

dsgelab --example rbc        # resolve um modelo do pacote e imprime os diagnósticos

Exemplos

Os exemplos funcionam sem alteração, na ordem em que aparecem, logo depois do pip install. Um verificador (tools/check_site_examples.py) extrai e executa cada bloco desta página, de modo que um bloco que deixe de rodar faz a verificação falhar. Artigos replicados no pacote, com o número que a página confere, estão em Replicações.

1. Resolver um modelo que já vem no pacote

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)     # desvios em %
mom = sol.moments(cov={"eps_a": 0.007})           # dp, correlações, autocorrelações
dec = sol.variance_decomposition(cov={"eps_a": 0.007})

2. Declarar o seu próprio modelo

O modelo é declarado com três listas de nomes e uma função, e as variáveis são acessadas por nome, em cur para o período t e nxt para 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),                    # recursos
        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,              # tecnologia
    ]

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

sol_meu.P                     # transição dos estados
sol_meu.accuracy()            # a ordem implícita tem de dar ~2

Quando as variáveis são lidas de um vetor pela posição, trocar um índice por outro não produz mensagem de erro, e o programa resolve, sem aviso, um modelo diferente do pretendido. Com cada variável lida pelo nome, essa troca deixa de ser possível.

3. As mesmas equações, em texto, com derivadas exatas

from dsgelab.symbolic import from_equations

simbolico = 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_sim = dsg.solve(simbolico, P(), guess={"k": 0.2, "a": 1.0, "c": 0.36})

dsg.compare_irfs(sol_meu, sol_sim, "eps_a")["worst"]   # os dois são o mesmo modelo

A jacobiana e a hessiana do resíduo saem exatas, por diferenciação simbólica. A solução de segunda ordem é construída sobre a hessiana, e uma diferença finita de segunda ordem divide pelo quadrado do passo, de modo que uma hessiana numérica conserva cerca de metade dos dígitos significativos do resíduo a partir do qual foi calculada.

4. Segunda ordem e correções de risco

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

so.g_ss, so.h_ss          # o que a primeira ordem põe em zero
so.mean_correction()      # média ergódica, em forma fechada
so.ergodic_mean()         # a mesma coisa, em nível
so.accuracy()             # a ordem implícita tem de dar ~3
so.simulate(periods=5000, seed=1)   # podada, por padrão

Uma simulação de primeira ordem mantém, por construção, a média dos estados no estado estacionário. Quando o deslocamento das regras de decisão é grande, ela pode errar a magnitude e também o sinal de um efeito sobre a média, e o ordenamento de políticas pelo bem-estar depende desse sinal.

5. Transição sob previsão perfeita, com política anunciada

from dataclasses import replace

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

T = 80
surpresa = dsg.perfect_foresight(meu, ss0, ss1, params=p_novo, T=T)

cronograma = [p_velho] * 10 + [p_novo] * (T + 1 - 10)      # anunciado em t=0
anunciada = dsg.perfect_foresight(meu, ss0, ss1, params=cronograma, T=T)

surpresa["k"], anunciada["k"]              # níveis, (T+1,)
anunciada["c"][0] != surpresa["c"][0]      # o anúncio aparece no impacto

A transição é resolvida pelo método de Newton com os períodos empilhados (stacked time), sobre a mesma função de resíduos que a perturbação lineariza. Uma surpresa entra como mudança em terminal com parâmetros constantes, e um cronograma legal conhecido de antemão entra como sequência de valores dos parâmetros, um por período.

6. Calibrar parâmetros por alvos e saber se estão identificados

def observaveis(ss, p):
    y = ss["a"] * ss["k"] ** p.alpha
    return {"ky": ss["k"] / y, "iy": p.delta * ss["k"] / y}   # capital-produto, participação do investimento

res = dsg.calibrate(
    meu, 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=observaveis,
    bounds={"beta": (0.5, 0.999), "delta": (0.001, 1.0)},
)
print(res.report())

ident = dsg.target_jacobian(meu, replace(P(), **res.values), free=["beta", "delta"],
                            targets=["ky", "iy"], guess={"k": 0.2, "a": 1.0, "c": 0.36},
                            observables=observaveis)
ident["identified"]        # True: dois alvos fixam os dois parâmetros

um = dsg.target_jacobian(meu, P(), free=["beta", "alpha"], targets=["k"],
                         guess={"k": 0.2, "a": 1.0, "c": 0.36})
um["identified"]           # False: dois parâmetros, um alvo
um["null_directions"]      # a direção que nenhum alvo enxerga

A razão capital-produto e a participação do investimento determinam β e δ em conjunto. Em forma fechada, δ = iy/ky = 0,025 e β = 1/(α/ky + 1 − δ) = 0,98912, que é o que calibrate devolve. Avaliada nessa solução, target_jacobian informa que os dois alvos identificam os dois parâmetros. Com β e α livres e o estoque de capital como único alvo, a matriz de elasticidades tem posto um, e null_directions dá a combinação dos dois parâmetros que deixa o alvo inalterado. Feita à mão, a calibração vira bisseções aninhadas, uma por parâmetro, e cada uma trata o seu alvo como se dependesse só do seu parâmetro, embora os alvos sejam determinados conjuntamente. Quando essa matriz não tem posto completo, a calibração não está identificada, e apertar a tolerância do solver não muda isso.

7. Verificações independentes

dsg.blanchard_kahn_report(sol)["margin"]        # a folga, não só o "passou"

dsg.check_slack(sol, {"c_positivo": lambda v: v["c"]}, size=-0.5)   # ao longo das IRFs

dsg.steady_state_invariance(                    # uma invariância que é teorema do modelo
    lambda v: dsg.steady_state(meu, {"k": 0.2, "a": 1.0, "c": 0.36}, replace(P(), rho=v)),
    [0.5, 0.8, 0.95], ["k", "c"])["ok"]         # k não depende da persistência

8. Derivar as condições de primeira ordem (CPOs) do problema declarado

Basta declarar o problema, e o pacote diferencia o lagrangiano. Como as derivadas são calculadas simbolicamente, as condições ficam livres dos sinais trocados e dos termos esquecidos da diferenciação à mão, erros que não alteram o estado estacionário e por isso passam despercebidos quando só ele é conferido.

from dsgelab.blocks import Problem, optimality_conditions

prob = Problem(
    "log(c)",                                   # retorno do período
    controls=["c", "i"], states=["k"], exogenous=["a"],
    constraints=["c + i = a*k**alpha",          # recursos
                 "k(+1) = (1-delta)*k + i"],    # acumulação
    discount="beta", parameters=["alpha", "beta", "delta", "rho"],
)
foc = optimality_conditions(prob)

for eq in foc["equations"]:
    print(f"{eq:60s}  # {foc['detail'][eq]}")   # cada equação diz de onde veio

derivado = 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"])

A equação de Euler de cada estado sai da derivada do lagrangiano do período seguinte em relação ao estado herdado. Condições de transversalidade, restrições de desigualdade e problemas cuja condição intertemporal não tem essa forma ficam a cargo do modelador. O pacote deriva as condições do problema declarado, de modo que um erro na declaração passa para as equações. Escolher o problema continua sendo tarefa de quem modela.

9. Terceira ordem

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

to.g_kkk, to.h_kkk        # os termos cúbicos
to.accuracy()             # a ordem implícita tem de dar ~4
to.simulate(periods=5000, seed=1)    # podada em três camadas

10. Exportar para o Dynare e para o MacroModelling.jl

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

texto = dsg.to_dynare(simbolico, P(), steady=ss_sim, cov={"eps_a": 0.007},
                     path="modelo.mod", order=2)
jl = dsg.to_macromodelling(simbolico, P(), cov={"eps_a": 0.007}, path="modelo.jl")

"predetermined_variables k;" in texto     # a convenção temporal é traduzida automaticamente

No dsgelab o estado em t é o disponível em t. No Dynare a variável predeterminada é a decidida em t. Uma conversão manual tem de traduzir uma convenção na outra. O exportador faz essa tradução e também escreve os processos exógenos um período atrás, para que a inovação apareça sem defasagem. Na ida e volta, o pacote escreve o .mod e o Dynare o resolve. Nos três modelos de referência, as duas soluções diferem em menos de 1e-13 no estado estacionário, nas IRFs e nos momentos.

11. Política ótima pelo planejador de Ramsey

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"],                          # choque de custo: processo exógeno
    states=["u"], jumps=["x", "pi", "i"], shocks=["e"],
    parameters=["beta", "sigma", "kappa", "lam", "rho"],
    instruments=["i"],                              # a variável sem equação própria
    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   # regra de metas

As condições do planejador são derivadas simbolicamente das equações do setor privado, e a função devolve um Model comum, com os multiplicadores e as suas defasagens como variáveis. A perspectiva é a atemporal, a mesma do ramsey_model do Dynare. No modelo novo-keynesiano a resposta segue a regra de metas de Galí, π_t = −(λ/κ)(x_t − x_{t−1}), com erro abaixo de 1e-12.

12. Política discricionária, regra simples ótima e bem-estar

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"]                       # sem compromisso
w_comp = dsg.welfare(nk, pn, steady=ss_nk, cov={"e": 0.01}).unconditional
w_disc < w_comp < 0                                            # o compromisso vale algo

regra = 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",           # a regra simples
                           "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(regra, 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

A política discricionária segue o algoritmo de Dennis (2007), o mesmo do discretionary_policy do Dynare, e reproduz a solução fechada de Galí para o modelo novo-keynesiano. A regra simples minimiza a perda quadrática nos momentos incondicionais, como o osr do Dynare. O bem-estar de segunda ordem, condicional e incondicional, é o do evaluate_planner_objective.

13. Choques antecipados e surpresas na previsão perfeita

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)

anunciado = np.zeros((80, 1)); anunciado[9, 0] = 0.02         # sabido desde o início
antecipada = dsg.perfect_foresight(rbc_pf, ss_pf, ss_pf, params=ppf, T=80, shocks=anunciado)
surpresa = dsg.perfect_foresight_with_expectation_errors(       # só se sabe em t = 10
    rbc_pf, ss_pf, ss_pf, params=ppf, T=80,
    news=[dsg.News(period=10, shocks=anunciado)])

antecipada["c"][0] != ss_pf.values["c"]      # o consumo reage ao anúncio
abs(surpresa["c"][0] - ss_pf.values["c"]) < 1e-12   # e não reage ao que ninguém sabe

Um choque antecipado (um objeto News) é uma nova trajetória esperada de choques e parâmetros, com um novo estado estacionário terminal, e a cada um os agentes replanejam a partir dos estados já realizados. O método é o do perfect_foresight_with_expectation_errors_solver do Dynare, e ModFile.news converte um .mod com a tabela de informação em CSV na lista de choques antecipados (objetos News).

14. Estimação bayesiana

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))
serie = sol_ar.to_levels(sol_ar.simulate(periods=400, cov={"e": 0.1}, seed=1))["y"]

pb = dsg.EstimationProblem(
    model=ar, data={"y": serie}, 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))
moda = dsg.posterior_mode(pb)                       # moda, hessiana, Laplace
cadeia = dsg.metropolis_hastings(pb, moda, draws=1500, chains=2, seed=1)

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

O EstimationProblem junta o modelo, os dados, as distribuições a priori e os observáveis. Essas distribuições seguem as do Dynare, definidas por média e desvio-padrão. Os observáveis podem estar em nível ou em log, ter erro de medida e ter dados faltantes. A partir do problema, o pacote calcula a moda a posteriori com a aproximação de Laplace, a máxima verossimilhança, o Metropolis com R-hat e tamanho efetivo da amostra, o SMC, a verificação preditiva a priori e a verossimilhança marginal. ModFile.estimation_problem constrói um EstimationProblem a partir de um .mod com estimated_params e varobs.

15. Restrições ocasionalmente ativas (OccBin)

credito = dsg.OccbinModel.from_equations(
    ["c = y + b - R*LB",
     "LB(+1) = b",
     "[name='divida', relax='folga'] b = m*y",      # o limite vale ...
     "[name='divida', bind='folga'] lam = 0",       # ... ou o multiplicador é zero
     "lam = 1/c - beta*R/c(+1)",
     "log(y(+1)) = rho*log(y) + u"],
    constraints={"folga": {"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(credito, 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["folga"].any()                                   # a renda alta afrouxa o limite
max(abs(ob.piecewise["lam"][ob.binding["folga"]])) < 1e-14  # e o multiplicador vai a zero

Cada regime é linearizado no estado estacionário de referência, com o seu resíduo como constante. Dadas as datas em que cada restrição vale, a trajetória segue regras que variam no tempo, calculadas de trás para a frente, e a conjectura inicial sobre essas datas é revista pela regra do +occbin do Dynare, com curb_retrench se pedido. ModFile.occbin_simulate() simula um .mod com equações marcadas (bind=, relax=), occbin_constraints e shocks(surprise).

16. Trajetória estendida, determinística e estocástica

arvore = dsg.from_equations(                       # Burnside (1998): preço de uma árvore de Lucas
    ["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(arvore, dict(x=0.0179, y=g0/(1 - g0)), pa)
sorteio = 0.0348*np.random.default_rng(0).standard_normal((20, 1))

ep = dsg.extended_path(arvore, ss_a, pa, sorteio, horizon=60)                 # sem risco à frente
sep = dsg.extended_path(arvore, ss_a, pa, sorteio, horizon=60, order=1,       # risco do próximo período
                       cov=[[0.0348**2]])
bool((sep["y"][1:] > ep["y"][1:]).all())       # com theta < 0, o risco encarece a árvore

A cada período chega o choque, calculam-se os estados a que ele leva e os agentes resolvem um problema não linear que alcança o estado estacionário horizon períodos depois. Com order > 0 os choques dos próximos períodos são integrados por quadratura de Gauss–Hermite numa árvore de planos, esparsa ou cheia, com a correção de risco de segunda ordem opcional (hybrid). O algoritmo e essas opções são os do extended_path do Dynare.

17. Declaração do modelo: defasagens, calibração no modelo, homotopia, diagnóstico

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)"],                         # a defasagem de um salto vira um estado auxiliar
    states=["k"], jumps=["c", "q", "dq"], shocks=[], parameters=["alpha", "beta", "delta"],
    calibration_equations={"beta": "k/q = 10"})  # beta é o que põe k/q em 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)})   # a conjectura inicial serve para delta = 0,02

Um estado estacionário em forma fechada (steady_state_model, aqui ou num .mod) só é usado depois de conferido com as equações, e é recusado, como no Dynare, se não as zera. dsg.assemble junta os problemas de vários agentes e as condições de mercado num só modelo.

18. Tendências e crescimento balanceado

tendencia = 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 = tendencia.steady_state()                  # o modelo já sai destendenciado

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

Num .mod escrito em níveis, as variáveis declaradas com trend_var e var(deflator=) são destendenciadas como no Dynare, e uma equação que continue dependendo do nível da tendência é recusada. balanced_growth calcula os fatores de crescimento e as razões entre variáveis que compartilham a tendência, o mesmo problema que o bgp.write do Dynare resolve.

O que o pacote oferece

A única convenção que você precisa acertar

eps é a inovação datada em t+1, desconhecida em t, com E_t eps = 0. Um AR(1) se escreve assim

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

Validação

A validação tem três camadas. O conjunto de testes compara o pacote com soluções fechadas e respostas conhecidas analiticamente. Os scripts de validation/ o comparam com o Dynare nos exemplos e arquivos de teste do próprio Dynare, com o MacroModelling.jl numa equação de calibração e com o código MATLAB original de Heutel (2012), apresentado em Replicações. Uma última comparação, num modelo aplicado de um artigo não publicado do autor, não pode ser refeita a partir do repositório.

Na segunda camada, tanto os modelos quanto os programas de referência são de outros autores. A terceira usa os mesmos programas de referência num modelo do próprio autor. Cada recurso é validado com uma referência independente quando ela existe (outra implementação, uma solução fechada ou um resultado publicado). Alguns se apoiam só em testes internos, a saber as autocorrelações, a decomposição da variância, a poda (pruning), a IRF generalizada de segunda ordem, as trajetórias de parâmetros anunciadas e a verificação de Blanchard–Kahn, e a exportação para o MacroModelling.jl não tem teste de ida e volta.

Comparação com soluções analíticas

Dos 148 testes do pacote, os desta camada comparam o resultado com respostas que podem ser deduzidas no papel. Com utilidade logarítmica e depreciação total, a solução do RBC tem de ser igual a P = [[α, 1], [0, ρ]]. No modelo novo-keynesiano de três equações com φx = 0,125, a regra de Taylor tem de dar equilíbrio determinado com φπ = 1,5 e indeterminado com φπ = 0,8. Num AR(1) lognormal, o coeficiente de segunda ordem h_ss tem de valer sigma², de modo que a regra de decisão se desloca em sigma²/2, o termo de Jensen. A transição sob previsão perfeita é comparada com o caminho fechado k' = αβk^α. O teste da ordem de acurácia exige que dividir o deslocamento pela metade divida o resíduo de primeira ordem por quatro e o de segunda por oito.

Comparação com os modelos de referência do Dynare

Os exemplos de dynare.org/resources foram executados no Dynare e comparados com o dsgelab variável a variável.

Modelo	Estado estacionário	IRFs	Momentos teóricos
example1 (Collard, choques correlacionados)	1,2e-14	2,9e-14	2,6e-14	OK
example2 (o mesmo em logaritmos)	2,1e-05	4,3e-06	3,3e-06	OK*
RBC_baseline (Pfeifer, 15 variáveis)	0,0e+00	2,3e-14	3,5e-14	OK

* No example2, o steady do Dynare encerra a busca cedo, num ponto que deixa resíduo de 1,05e-06 nas equações do modelo, contra 1,03e-44 do dsgelab, e toda a discrepância decorre desse resíduo. Como example1 e example2 são a mesma economia em coordenadas diferentes, eles permitem uma verificação independente, e exp(e.e. do example2) coincide com o example1 em nível até 1e-15.

Os exemplos do Dynare, lidos diretamente, sem conversão manual

dsg.from_mod lê o .mod como está, com macros, variáveis locais, avanços e defasagens de qualquer ordem e choques dentro de qualquer equação. O estado estacionário, as IRFs e os desvios-padrão teóricos que saem dele são comparados com os que o Dynare calcula a partir do mesmo arquivo.

Modelo	Variáveis	Estado estacionário	IRFs	Momentos
fs2000 (Schorfheide 2000)	14	9,2e-15	3,4e-13	3,4e-13	OK
NK_baseline (Fernández-Villaverde e Rubio-Ramírez)	28	9,3e-14	5,9e-11	2,2e-11	OK
Gali_2015 (cap. 3, raiz unitária)	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 (em níveis, trend_var e deflator)	11	1,4e-08	1,5e-08	2,2e-08	OK*

* A concordância é limitada pelo resíduo de 4e-07 que o estado estacionário do Dynare deixa nas equações (7e-08 no fs2000_nonstationary, que cada um dos dois programas destendencia por si). Os três modelos da tabela anterior também passam por este caminho. No example1 e no RBC_baseline os desvios ficam abaixo de 1,2e-13, e o maior é o do estado estacionário do RBC_baseline, 1,1e-13. O example2 dá os números daquela tabela.

Política ótima comparada com o ramsey_model do Dynare

A comparação cobre o estado estacionário e as IRFs das variáveis privadas e dos multiplicadores, que os dois programas numeram pela mesma equação e com o mesmo sinal.

Modelo	Estado estacionário	IRFs
nk_linear (Galí, cap. 5)	0,0e+00	5,3e-15	OK
rbc_labortax (capital como estoque escolhido)	5,0e-11	4,9e-13	OK
nk_ramsey_model (testes do Dynare, não linear)	2,7e-12	2,1e-14	OK

Política discricionária, regras simples e bem-estar comparados com o Dynare

from_mod lê cada arquivo sem alteração. Três vêm dos testes do Dynare, e nk_osr_planner.mod e rbc_labortax.mod foram escritos neste repositório.

Teste	O que se compara	Desvio
osr_example.mod	perda no ótimo do Dynare	2,6e-16	OK*
nk_osr_planner.mod	regra simples ótima sobre o bem-estar	1,5e-16	OK
dennis_1.mod	IRFs sob política discricionária; bem-estar incondicional e condicional	4,3e-08; 5,3e-10, 6,1e-10	OK
nk_ramsey_model.mod	Ramsey lido do .mod; bem-estar de 2ª ordem	2,6e-14; 1,4e-12	OK
rbc_labortax.mod	bem-estar de 2ª ordem com capital	3,5e-13	OK

* O ótimo do dsgelab tem perda de 0,4266626, menor que os 0,4268336 do Dynare. Avaliada pelo próprio Dynare no ponto do dsgelab, a perda difere da do dsgelab em 4,3e-15, e portanto a diferença vem do otimizador do Dynare, que encerra a busca cedo numa superfície plana. No ótimo, a regra neutraliza o choque de demanda, gamma_y_ = 1/|σ| = 16,667. O piso de 4e-08 nas IRFs sob política discricionária é a tolerância do Dynare, sqrt(eps).

Previsão perfeita comparada com o Dynare

Os arquivos de tests/deterministic_simulations do Dynare são lidos por from_mod, e a comparação cobre toda variável em todo período.

Teste	O que testa	Desvio
rbc_det1.mod	condição inicial (histval)	1,1e-14	OK
rbc_det2.mod	choque no primeiro período	6,6e-15	OK
rbc_det3.mod	choques antecipados, períodos 4 a 8	8,0e-14	OK
pfwee.mod	erros de expectativa, choques transitórios e permanentes	2,8e-12	OK
pfwee_constant_sim_length.mod	idem, horizonte constante	3,0e-12	OK

Por padrão, o Dynare interrompe a busca quando o resíduo chega a 1e-5 e recalcula o estado estacionário terminal com tolerância eps^(1/3). Com esses valores padrão, o k terminal do pfwee difere da solução fechada em 4e-5. A validação executa o Dynare com 1e-12.

Estimação comparada com o Dynare

O fs2000.mod (Schorfheide, 2000) da pasta examples/ do Dynare é lido por from_mod com as suas distribuições a priori e os seus dados.

O que se compara	Desvio
verossimilhança e densidade a priori, no ponto inicial	6e-12; 4e-13	OK
variáveis suavizadas; choques suavizados	8,4e-15; 1,1e-15	OK
shock_decomposition e realtime_shock_decomposition	5e-15	OK
previsão de 8 períodos; bandas de 90%	2e-15; 6e-09	OK
moda a posteriori (núcleo do dsgelab ≥ o do Dynare)	+6e-06	OK
Laplace com a covariância do Dynare	2e-13	OK
médias do Metropolis, em erros de Monte Carlo (o maior é 2,89)	< 3	OK
filtro difuso (nível de preços observado), em dois pontos	2,5e-12	OK
previsão condicional (tests/conditional_forecasts)	1,6e-13	OK

As verossimilhanças coincidem até a última casa decimal só com o riccati_tol do Dynare desligado. Na configuração padrão, o Dynare passa a usar o ganho estacionário do filtro quando a variância converge a 1e-6, o que altera o número na sétima casa decimal. A aproximação de Laplace calculada com a hessiana do próprio dsgelab difere da do Dynare em 0,03, diferença que vem das duas hessianas numéricas.

Não linearidades e declaração do modelo comparadas com o Dynare

from_mod lê os arquivos de tests/occbin, tests/ep, tests/bgp e tests/trend_var do Dynare 7.2, e cada variável é conferida em cada período.

Teste	O que se compara	Desvio
borrcon.mod	OccBin, uma restrição, curb_retrench: trajetórias por partes e linear; regimes	2,2e-16; 30/30	OK
dynrbc.mod	OccBin, duas restrições e quatro regimes, 52 choques não antecipados	1,6e-15; 169/169	OK
burnside.mod	trajetória estendida; estocástica de ordem 1 e 2, árvores esparsa e cheia, híbrida	1,1e-11	OK*
linearmodel0.mod	trajetória estendida num modelo linear; solução de primeira ordem do Dynare	2,0e-15; 3,1e-14	OK
solow, ramsey, nk	o problema do bgp.write: resíduos em cinco pontos; asserções dos testes	2,2e-16; 2,8e-12	OK
walras.mod, bad_steady.mod, homotopy.mod	model_diagnostics; steady_state_model recusado; homotopia	mesmo relatório	OK

* Aqui o Dynare roda com tolerância do solver de 1e-12, contra o valor padrão de 1e-5, e é essa tolerância que limita a concordância.

Uma equação de calibração comparada com o MacroModelling.jl

Num pequeno modelo RBC, a meta de estado estacionário k/q = 2.5 determina β, e o modelo foi resolvido pelo dsgelab e pelo MacroModelling.jl. Os dois devolvem β = 0,847457627118644, o valor da solução fechada, e os valores de estado estacionário diferem em no máximo 4,4e-15.

Um modelo aplicado do autor, comparado com o Dynare e o MacroModelling.jl

Num DSGE de 71 variáveis, três grupos de renda e cinco choques, as correções de média de segunda ordem do dsgelab diferem das do Dynare em até 1,1e-05 ponto percentual e das do MacroModelling.jl em até 1,2e-05, na calibração de outubro de 2026. O Dynare usa segunda ordem não podada e o MacroModelling, podada, enquanto o dsgelab resolve a equação de Sylvester e devolve a média em forma fechada, sem simular. O modelo é de um artigo não publicado e não está neste repositório, de modo que esta comparação, ao contrário das anteriores, não pode ser refeita a partir daqui.

Comparação com outros pacotes

Outros pacotes em Python resolvem modelos DSGE, e a tabela põe as capacidades de cada um ao lado das do dsgelab.

Capacidade	dsgelab	gEconpy	dolo	PyMacLab	econpizza
Deriva as CPOs do problema do agente	sim	sim	—	—	—
Perturbação de 1ª ordem	sim	sim	sim	sim	—
2ª ordem + correções de risco	sim	não	sim	sim	—
3ª ordem (bloco determinístico)	sim	—	sim	—	—
3ª ordem com risco na inclinação	sim (g_ssk, h_ssk), sem o tensor completo do Dynare	—	sim	—	—
Métodos globais / projeção	parcial (colocação de Chebyshev, até três estados)	—	sim	—	parcial
Agentes heterogêneos	—	—	—	—	sim
Previsão perfeita não linear	sim	sim	sim	?	sim
Estimação bayesiana	sim (moda, Metropolis–Hastings e SMC, com o filtro de Kalman sobre a solução de primeira ordem)	sim	—	—	—
Calibração por alvos + identificação	sim	parcial	?	?	?
Exportar para Dynare e MacroModelling.jl	sim	só Dynare	só Dynare	—	—
Ordem de acurácia autoverificada	sim (ordens 2, 3, 4)	?	?	?	?
Dependências	numpy, scipy	sympy, PyMC	pandas (versão antiga fixada)	numpy, scipy	JAX
Instala no Python 3.14	sim	sim	não	não	?

Levantamento de 9 de setembro de 2026, revisado em 11 e 12 de outubro de 2026, com o gEconpy 2.2. Versões posteriores desses pacotes podem mudar algumas células, e o docs/comparison.md indica a fonte de cada uma.

No dsgelab, um problema declarado vai das condições de primeira ordem derivadas a uma solução e a uma exportação para o Dynare ou o MacroModelling.jl, com a convenção temporal traduzida automaticamente. A mesma função de resíduos serve à perturbação e à transição determinística. A solução de segunda ordem é calculada dentro do próprio Python, sem chamar o Dynare pelo Octave. A condição de Blanchard–Kahn e a convergência do estado estacionário são verificadas a cada solução, a ordem de acurácia e a identificação de uma calibração exigem uma chamada cada, e as únicas dependências são numpy e scipy. Nos três usos citados na nota acima (estimação com o PyMC, métodos globais e agentes heterogêneos), os pacotes ali indicados oferecem mais que o dsgelab. A tabela completa, que separa as células verificadas com a instalação de um pacote das que foram lidas só na documentação, está em docs/comparison.md.

Como citar

Se o pacote contribuiu para um resultado publicado, cite-o nas referências como qualquer outro insumo metodológico, e informe a versão, porque um número que só se reproduz numa versão específica não é reproduzível sem ela.

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    = {Pacote em Python}
}

ABNT

MOREIRA, Gustavo Hoffmann. dsgelab: a general toolkit for solving DSGE models in Python.
Versão 0.5.0. 2026. Software em Python. DOI: 10.5281/zenodo.23456475. Disponível em:
https://doi.org/10.5281/zenodo.23456475. Acesso em: [data].

No corpo do texto

Cite o nome e a versão junto do resultado, como no exemplo abaixo.

O modelo foi resolvido com o dsgelab 0.5.0 (Moreira, 2026), em perturbação de segunda ordem.

Como citar as ferramentas usadas pelo pacote

Se o seu resultado dependeu delas, cite-as também. Cite o Dynare quando o modelo foi exportado para ele ou conferido nele, e o MacroModelling.jl na mesma situação. NumPy e SciPy pedem citação explicitamente, e todo cálculo do dsgelab roda sobre elas.

Documentação