DSGE Estimation

MacroEconometricModels.jl provides two paradigms for estimating the deep structural parameters of DSGE models. Frequentist estimation via estimate_dsge matches model-implied moments to data moments using Generalized Method of Moments (GMM) with four moment conditions. Bayesian estimation via estimate_dsge_bayes combines prior distributions with the likelihood function, targeting the posterior with Sequential Monte Carlo (SMC), SMC$^2$, or Random-Walk Metropolis-Hastings (RWMH). Both approaches build on the solution infrastructure documented in DSGE Models.

using MacroEconometricModels, Random, Distributions
Random.seed!(42)
# Simple DSGE for fast estimation
spec = @dsge begin
    parameters: rho = 0.9, phi = 1.5
    endogenous: y, i
    exogenous: e
    y[t] = rho * y[t-1] + e[t]
    i[t] = phi * y[t]
end
spec = compute_steady_state(spec)
sol = solve(spec; method=:gensys)
Y_data = simulate(sol, 200)

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

Quick Start

Recipe 1: IRF matching GMM

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

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

Recipe 2: Bayesian SMC

result = estimate_dsge_bayes(spec, Y_data, [0.9];
    priors=Dict(:rho => Beta(5, 2)),
    method=:smc, observables=[:y], n_smc=50)
report(result)

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

Recipe 3: SMC$^2$ with projection solver

result = estimate_dsge_bayes(spec, Y_data, [0.9];
    priors=Dict(:rho => Beta(5, 2)),
    method=:smc2, observables=[:y], n_smc=200, n_particles=100,
    solver=:projection, solver_kwargs=(degree=5,))
report(result)

Recipe 4: Bayesian IRFs with credible bands

# Dual 68%/90% credible bands from posterior draws
birf = irf(result, 20; n_draws=10)
report(birf)

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

plot_result(birf)

Recipe 5: Bayesian FEVD with credible bands

bfevd = fevd(result, 20; n_draws=10)
report(bfevd)

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

Recipe 6: Model comparison via Bayes factor

log_bf = bayes_factor(result1, result2)

A positive value favors model 1; following Kass & Raftery (1995), $2 \cdot \log \text{BF} > 6$ constitutes strong evidence.

GMM Estimation

estimate_dsge estimates deep structural parameters by matching model-implied moments to data moments via Generalized Method of Moments. Four moment conditions are available: IRF matching, Euler equation GMM, simulated method of moments (SMM), and analytical GMM.

IRF Matching (Christiano, Eichenbaum & Evans 2005)

The IRF matching estimator minimizes the distance between model-implied and empirical impulse response functions:

\[\hat{\theta} = \arg\min_\theta \; \big(\Phi^m(\theta) - \Phi^d\big)' \, W \, \big(\Phi^m(\theta) - \Phi^d\big)\]

where:

• $\Phi^m(\theta)$ is the vector of stacked model-implied IRFs at parameter $\theta$
• $\Phi^d$ is the vector of empirical IRFs estimated from a VAR(p) on the data
• $W$ is the GMM weighting matrix (identity for one-step, inverse of $\hat{\Omega}$ for two-step)

The procedure first estimates a reduced-form VAR on the observed data, computes Cholesky-identified IRFs, then searches over the structural parameter space to find the $\theta$ that best replicates those empirical IRFs. This is the workhorse approach for medium-scale DSGE estimation in the frequency domain.

est_irf = estimate_dsge(spec, Y_data, [:rho];
                    method=:irf_matching, var_lags=4, irf_horizon=20,
                    weighting=:two_step)
report(est_irf)

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

The two-step weighting uses the inverse of the estimated variance of the moment conditions from the first-step residuals. For pre-computed target IRFs (e.g., from a sign-identified VAR), pass them via target_irfs to bypass the internal VAR estimation.

Euler Equation GMM (Hansen & Singleton 1982)

The Euler equation approach exploits the model's optimality conditions directly as moment conditions:

\[E\Big[f\big(y_t, y_{t-1}, y_{t+1}, \varepsilon_t, \theta\big) \otimes z_t\Big] = 0\]

where:

• $f(\cdot)$ is the vector of Euler equation residuals from the DSGE specification
• $z_t$ is a vector of instruments (lagged endogenous variables)
• $\otimes$ denotes the Kronecker product forming the interaction of residuals and instruments

This method does not require solving the model –- it evaluates the equilibrium conditions directly on the data. The instrument set consists of $n_{\text{lags}}$ lags of the endogenous variables, producing $n_{\text{eq}} \times n_{\text{vars}} \times n_{\text{lags}}$ moment conditions.

est_euler = estimate_dsge(spec, Y_data, [:rho];
                    method=:euler_gmm, n_lags_instruments=4,
                    weighting=:two_step)
report(est_euler)

Simulated Method of Moments

When analytical moments are unavailable, SMM matches sample moments from simulated data to their empirical counterparts. The simulation ratio $S/T$ (default: 5) controls the trade-off between computational cost and simulation noise:

\[\hat{\theta} = \arg\min_\theta \; \big(\hat{m}_S(\theta) - \hat{m}_T\big)' \, W \, \big(\hat{m}_S(\theta) - \hat{m}_T\big)\]

where:

• $\hat{m}_S(\theta)$ is the vector of moments computed from a simulated path of length $S$
• $\hat{m}_T$ is the vector of sample moments from the observed data of length $T$
• $W$ accounts for both sampling and simulation uncertainty

est_smm = estimate_dsge(spec, Y_data, [:rho];
                    method=:smm, sim_ratio=5)
report(est_smm)

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

The default moment function computes autocovariances at lag 1. Supply a custom moments_fn to target specific features of the data.

Analytical GMM

Analytical GMM computes model-implied moments from the unconditional distribution without simulation. For linear models, the Lyapunov equation provides exact second moments via analytical_moments. For higher-order perturbation solutions, moments are computed from pruned simulations.

est_agmm = estimate_dsge(spec, Y_data, [:rho];
                    method=:analytical_gmm, solve_method=:gensys,
                    solve_order=1, lags=1)
report(est_agmm)

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

Use solve_method=:perturbation with solve_order=2 to match moments from second-order solutions, which capture risk premia and precautionary behavior absent from linear approximations.

Hansen J-test

When the number of moment conditions exceeds the number of parameters, the Hansen (1982) J-statistic tests whether the over-identifying restrictions hold:

\[J = T \cdot g(\hat{\theta})' \, \hat{W} \, g(\hat{\theta}) \sim \chi^2(q - p)\]

where:

• $q$ is the number of moment conditions
• $p$ is the number of estimated parameters
• $g(\hat{\theta})$ is the sample moment vector evaluated at the estimated parameters

A large J-statistic (low p-value) indicates model misspecification –- the model cannot simultaneously satisfy all moment conditions.

est.J_stat     # Hansen J-statistic
est.J_pvalue   # p-value under chi-squared distribution

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

GMM Keywords

GMM Return Value (DSGEEstimation{T})

StatsAPI interface: coef(est) returns the parameter vector, vcov(est) the covariance matrix, stderror(est) the standard errors, and dof(est) the number of estimated parameters.

Bayesian Estimation

Bayesian estimation combines prior distributions $\pi(\theta)$ with the data likelihood $\mathcal{L}(Y|\theta)$ to characterize the posterior distribution:

\[p(\theta | Y) \propto \mathcal{L}(Y | \theta) \cdot \pi(\theta)\]

where:

• $p(\theta | Y)$ is the posterior distribution over structural parameters
• $\mathcal{L}(Y | \theta)$ is the likelihood function (Kalman filter for linear models, particle filter for nonlinear models)
• $\pi(\theta)$ is the prior distribution encoding economic beliefs

Three sampling algorithms target the posterior: Sequential Monte Carlo (SMC), SMC$^2$, and Random-Walk Metropolis-Hastings (RWMH). SMC is the recommended default –- it produces the marginal likelihood as a by-product, handles multimodal posteriors, and parallelizes naturally.

Prior Specification

Priors are specified as a Dict{Symbol, Distribution} mapping parameter names to distributions from Distributions.jl. Parameter bounds are inferred automatically from the distribution support:

priors = Dict(
    :rho => Beta(5, 2)          # persistence: mean ≈ 0.71, support [0,1]
)
nothing # hide

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

Sequential Monte Carlo (Herbst & Schorfheide 2014)

SMC draws from a sequence of tempered distributions that bridge the prior to the posterior:

\[p_\phi(\theta) \propto \mathcal{L}(Y|\theta)^\phi \; \pi(\theta), \qquad \phi \in [0, 1]\]

where:

• At $\phi = 0$, particles are distributed according to the prior $\pi(\theta)$
• At $\phi = 1$, particles approximate the full posterior $p(\theta|Y)$
• The tempering schedule $0 = \phi_0 < \phi_1 < \cdots < \phi_S = 1$ bridges between the two

The algorithm proceeds in six steps:

1. Initialize: Draw $N$ particles from the prior: $\theta^{(i)} \sim \pi(\theta)$
2. Temper: Set the adaptive schedule $0 = \phi_0 < \phi_1 < \cdots < \phi_S = 1$ targeting a fixed ESS fraction
3. Reweight: At stage $s$, compute incremental weights $w^{(i)} \propto \mathcal{L}(Y|\theta^{(i)})^{\phi_s - \phi_{s-1}}$
4. Resample: If ESS falls below the threshold, resample particles via systematic resampling
5. Mutate: Apply $n_{\text{mh}}$ Metropolis-Hastings steps with proposal $q(\theta^*|\theta) = \mathcal{N}(\theta, \hat{\Sigma})$
6. Marginal likelihood: Estimate the normalizing constant: $\hat{p}(Y) = \prod_{s=1}^{S} \frac{1}{N} \sum_{i=1}^{N} w_s^{(i)}$

The adaptive tempering schedule selects $\phi_s$ to maintain the effective sample size at the target fraction (default: 50%). This avoids both degenerate weights (too large a step) and unnecessary computation (too small a step).

result_smc = estimate_dsge_bayes(spec, Y_data, [0.9];
    priors=Dict(:rho => Beta(5, 2)),
    method=:smc, observables=[:y], n_smc=50)
report(result_smc)

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

The likelihood is evaluated via the Kalman filter, which is exact for linear state-space models produced by Linear Solvers (:gensys, :blanchard_kahn, :klein).

SMC$^2$ (Chopin, Jacob & Papaspiliopoulos 2013)

SMC$^2$ nests a particle filter inside the outer SMC loop. At each mutation step, the likelihood $\mathcal{L}(Y|\theta^*)$ is evaluated by running a bootstrap particle filter rather than the Kalman filter. This enables Bayesian estimation of nonlinear DSGE models solved with Nonlinear Methods –- perturbation (order $\geq 2$), Chebyshev projection, or policy function iteration –- where the Kalman filter approximation breaks down.

The particle filter estimates the likelihood as:

\[\hat{\mathcal{L}}(Y|\theta) = \prod_{t=1}^{T} \frac{1}{M} \sum_{j=1}^{M} w_t^{(j)}\]

where $M$ is the number of inner particles (set via n_particles) and $w_t^{(j)}$ are the importance weights at time $t$.

result = estimate_dsge_bayes(spec, Y_data, [0.9];
    priors=Dict(:rho => Beta(5, 2)),
    method=:smc2, observables=[:y],
    n_smc=200, n_particles=100,
    solver=:projection, solver_kwargs=(degree=5,))
report(result)

The mutation step uses Conditional SMC (CSMC) to update both the parameter $\theta$ and the latent state trajectory jointly, maintaining a valid reference trajectory that prevents particle degeneracy.

Delayed Acceptance (Christen & Fox 2005)

Two-stage delayed acceptance pre-screens Metropolis-Hastings proposals with a cheap particle filter before running the expensive Conditional SMC. This preserves detailed balance while avoiding wasted computation on proposals that would be rejected.

Stage 1 (screening): Accept the proposal $\theta^*$ with probability

\[\alpha_1 = \min\!\Big(1, \; \exp\big[\phi \cdot \hat{\ell}_{\text{screen}}(\theta^*) + \log\pi(\theta^*) - \phi \cdot \hat{\ell}_{\text{screen}}(\theta) - \log\pi(\theta)\big]\Big)\]

Stage 2 (correction, only if Stage 1 accepts): Accept with probability

\[\alpha_2 = \min\!\Big(1, \; \exp\big[\phi \cdot \big(\hat{\ell}_{\text{full}}(\theta^*) - \hat{\ell}_{\text{screen}}(\theta^*)\big) - \phi \cdot \big(\hat{\ell}_{\text{full}}(\theta) - \hat{\ell}_{\text{screen}}(\theta)\big)\big]\Big)\]

where:

• $\hat{\ell}_{\text{screen}}$ is the log-likelihood from a bootstrap PF with $n_{\text{screen}}$ particles (cheap)
• $\hat{\ell}_{\text{full}}$ is the log-likelihood from CSMC with $n_{\text{particles}}$ particles (expensive)
• $\phi$ is the current tempering parameter

The product $\alpha_1 \cdot \alpha_2$ equals the standard MH acceptance probability in expectation, so the chain targets the exact posterior. The computational savings come from rejecting bad proposals cheaply at Stage 1 without ever running the full CSMC.

result = estimate_dsge_bayes(spec, Y_data, [0.9];
    priors=Dict(:rho => Beta(5, 2)),
    method=:smc2, observables=[:y],
    n_smc=200, n_particles=500,
    solver=:projection, solver_kwargs=(degree=5,),
    delayed_acceptance=true, n_screen=200)

Random-Walk Metropolis-Hastings

Standard Random-Walk Metropolis-Hastings with adaptive proposal scaling targeting the optimal 23.4% acceptance rate (Roberts & Rosenthal 2001). The proposal is a multivariate normal centered at the current draw:

\[\theta^* \sim \mathcal{N}\!\big(\theta^{(s)}, \; c \cdot \hat{\Sigma}\big)\]

where:

• $\hat{\Sigma}$ is the estimated posterior covariance (initialized from the prior, updated during burnin)
• $c$ is the step-size scalar adapted to target 23.4% acceptance

result_mh = estimate_dsge_bayes(spec, Y_data, [0.9];
    priors=Dict(:rho => Beta(5, 2)),
    method=:mh, observables=[:y],
    n_draws=50, burnin=25)
report(result_mh)

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

RWMH is simple to implement and diagnose but converges slowly for high-dimensional parameter spaces. For models with more than 5–10 parameters, SMC is strongly preferred.

Bayesian Keywords

Posterior Analysis

After estimation, five functions extract information from the BayesianDSGE result:

# Posterior summary: mean, median, std, 95% credible interval per parameter
ps = posterior_summary(result_smc)
ps[:rho][:mean]       # posterior mean of rho

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

ps[:rho][:ci_lower]   # lower bound of 95% CI for rho

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

# Log marginal likelihood (model evidence)
ml = marginal_likelihood(result_smc)

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

# Bayes factor: log p(Y|M₁) - log p(Y|M₂)
log_bf = bayes_factor(result1, result2)

# Prior vs posterior comparison table
tbl = prior_posterior_table(result_smc)
nothing # hide

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

# Posterior predictive simulation
Y_pred = posterior_predictive(result_smc, 10; T_periods=50)
size(Y_pred)

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

posterior_summary returns a Dict{Symbol, Dict{Symbol, T}} with keys :mean, :median, :std, :ci_lower (2.5th percentile), and :ci_upper (97.5th percentile) for each parameter. prior_posterior_table returns a vector of named tuples suitable for tabular display, comparing prior and posterior moments side by side. posterior_predictive draws n_sim parameter vectors from the posterior, solves the model at each, and simulates forward, returning an n_sim x T_periods x n_vars array of simulated paths.

Posterior IRFs and FEVD (Herbst & Schorfheide 2015)

Bayesian DSGE estimation quantifies parameter uncertainty. irf and fevd propagate this uncertainty into impulse responses and variance decompositions by re-solving the model at posterior parameter draws and computing pointwise credible bands.

For each of n_draws randomly selected posterior draws, the model is re-solved at those parameter values and the analytical IRF (or FEVD) is computed. The results are stacked and summarized with pointwise quantile bands. The default quantiles $[0.05, 0.16, 0.84, 0.95]$ produce dual 68% and 90% credible bands –- the standard reporting convention in the Bayesian DSGE literature.

# Bayesian IRFs with dual credible bands
birf_smc = irf(result_smc, 20; n_draws=10)
report(birf_smc)

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

plot_result(birf_smc)

# Bayesian FEVD
bfevd_smc = fevd(result_smc, 20; n_draws=10)
report(bfevd_smc)

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

# Custom quantiles (90% band only)
birf_90 = irf(result_smc, 20; n_draws=10, quantiles=[0.05, 0.95])
nothing # hide

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

Both methods return BayesianImpulseResponse{T} and BayesianFEVD{T} respectively –- the same types used by Bayesian VAR, so all existing report(), plot_result(), table(), and cumulative_irf() infrastructure works automatically.

Draws that produce indeterminate or explosive solutions are silently skipped. If all draws fail, an error is raised.

Posterior Predictive Simulation

simulate draws from the posterior predictive distribution with credible bands. For each posterior parameter draw, the model is re-solved and simulated forward T_periods periods:

bsim = simulate(result_smc, 50; n_draws=10)
report(bsim)

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

The result is a BayesianDSGESimulation{T} containing the pointwise median, quantile bands, and all raw simulation paths.

BayesianDSGESimulation{T} Return Value

Bayesian Return Value (BayesianDSGE{T})

StatsAPI interface: coef(result) returns the posterior mean parameter vector.

Complete Example

This example estimates the RBC model's persistence and shock volatility using both GMM and Bayesian methods, then compares the results:

using MacroEconometricModels, Random, Distributions
Random.seed!(42)

# 1. Specify the 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

# 2. Simulate data from the true model
sol = solve(spec)
Y_data = simulate(sol, 200)

# 3. Frequentist: IRF matching GMM
est_gmm = estimate_dsge(spec, Y_data, [:ρ, :σ];
                         method=:irf_matching, var_lags=4, irf_horizon=20)
report(est_gmm)

# 4. Bayesian: SMC estimation
priors = Dict(:ρ => Beta(5, 2), :σ => InverseGamma(2, 0.01))
est_bayes = estimate_dsge_bayes(spec, Y_data, [0.9, 0.01];
    priors=priors, method=:smc, observables=[:Y, :C], n_smc=2000)
report(est_bayes)

# 5. Posterior analysis
ps = posterior_summary(est_bayes)
tbl = prior_posterior_table(est_bayes)
ml = marginal_likelihood(est_bayes)

# 6. Bayesian IRFs with credible bands
birf = irf(est_bayes, 20; n_draws=100)
report(birf)

# 7. Bayesian FEVD
bfevd = fevd(est_bayes, 20; n_draws=100)
report(bfevd)

# 8. Posterior predictive simulation
bsim = simulate(est_bayes, 200; n_draws=100)
report(bsim)

The GMM point estimates and Bayesian posterior means should cluster near the true values $\rho = 0.9$ and $\sigma = 0.01$. The Bayesian IRFs propagate parameter uncertainty into impulse response bands, and the marginal likelihood enables formal model comparison.

Common Pitfalls

1. Wrong steady state: If the steady state is incorrect, all estimation methods fail silently –- the model solves to a nonsensical equilibrium, and the optimizer converges to economically meaningless parameter values. Always verify compute_steady_state and check that the solution satisfies is_determined(sol) before estimation.

2. Indeterminate model at prior draws: SMC initializes particles from the prior. If many prior draws produce indeterminate or explosive models, the likelihood evaluates to $-\infty$ and particles are wasted. Tighten priors to concentrate mass on the determinacy region, or increase n_smc to compensate.

3. Too few SMC particles: For posteriors with ridges, multimodality, or strong correlations, $n_{\text{smc}} = 1000$ may not suffice. Start with 5000+ and reduce only after confirming that the ESS history remains above the target and that repeated runs produce consistent marginal likelihood estimates.

4. Observable mismatch: The observables keyword specifies which endogenous variables in the model correspond to columns in the data matrix, in order. Mismatched dimensions or incorrect ordering produce nonsensical likelihood values. The number of observables must equal the number of data columns.

5. Solver choice for Bayesian estimation: Use :smc with :gensys (or :blanchard_kahn, :klein) for linear models –- the Kalman filter provides the exact likelihood. Use :smc2 with :projection or :pfi for nonlinear models where the particle filter is necessary. Using :smc with a nonlinear solver silently falls back to a first-order Kalman approximation that ignores higher-order dynamics.

References

• An, S., & Schorfheide, F. (2007). Bayesian Analysis of DSGE Models. Econometric Reviews, 26(2-4), 113-172. DOI

• Chopin, N., Jacob, P. E., & Papaspiliopoulos, O. (2013). SMC$^2$: An Efficient Algorithm for Sequential Analysis of State Space Models. Journal of the Royal Statistical Society: Series B, 75(3), 397-426. DOI

• Christen, J. A., & Fox, C. (2005). Markov Chain Monte Carlo Using an Approximation. Journal of Computational and Graphical Statistics, 14(4), 795-810. DOI

• Christiano, L. J., Eichenbaum, M., & Evans, C. L. (2005). Nominal Rigidities and the Dynamic Effects of a Shock to Monetary Policy. Journal of Political Economy, 113(1), 1-45. DOI

• Hansen, L. P. (1982). Large Sample Properties of Generalized Method of Moments Estimators. Econometrica, 50(4), 1029-1054. DOI

• Hansen, L. P., & Singleton, K. J. (1982). Generalized Instrumental Variables Estimation of Nonlinear Rational Expectations Models. Econometrica, 50(5), 1269-1286. DOI

• Herbst, E., & Schorfheide, F. (2014). Sequential Monte Carlo Sampling for DSGE Models. Journal of Applied Econometrics, 29(7), 1073-1098. DOI

• Kass, R. E., & Raftery, A. E. (1995). Bayes Factors. Journal of the American Statistical Association, 90(430), 773-795. DOI

• Roberts, G. O., & Rosenthal, J. S. (2001). Optimal Scaling for Various Metropolis-Hastings Algorithms. Statistical Science, 16(4), 351-367. DOI