Dynamic Stochastic General Equilibrium Models

Specify, solve, and estimate dynamic stochastic general equilibrium (DSGE) models by maximum likelihood and Bayesian methods. Supports both linear models via an equation-based formula interface and nonlinear models via string-based equations with perturbation up to third order (Schmitt-Grohe and Uribe, 2004 ). Solution uses the method of undetermined coefficients (Klein, 2000 ). Likelihood evaluated via the Kalman filter or a bootstrap particle filter (Gordon et al., 1993). Bayesian estimation uses adaptive Random-Walk Metropolis-Hastings or Particle Marginal Metropolis-Hastings (Andrieu et al., 2010 ) with parallel chain support. Additional tools include Bayes factor model comparison with Kass-Raftery evidence scales, Ramsey optimal policy via linear-quadratic regulator, nonlinear perfect foresight via stacked-time Newton (Juillard et al., 1998), Kalman smoothing, historical shock decomposition, local identification diagnostics, parameter sensitivity analysis, occasionally binding constraints, impulse-response functions, forecasting, and robust standard errors.

dsge

Dynamic Stochastic General Equilibrium Models for R.

Overview

The dsge package provides a comprehensive framework for specifying, solving, and estimating DSGE models entirely in R. No external software (Dynare, MATLAB, Octave) is required.

Key capabilities:

• Dynare model import (read_dynare()) -- read a Dynare .mod file straight into dsge (macro directives, calibration, steady state, shocks, measurement errors, priors, Ramsey/discretionary/OSR policy and OccBin constraints), and solve or estimate it without Dynare or MATLAB; MATLAB code in the file and _steadystate.m files are run by a built-in MATLAB interpreter, and perfect-foresight simulations (initval/endval, deterministic shocks, lmmcp constraints) by simulate_perfect_foresight()
• Linear models via formula interface (obs(), unobs(), state())
• Nonlinear models via string-based equations with perturbation up to third order (dsgenl_model(), solve_dsge(order = 1, 2, 3))
• Maximum likelihood estimation via the Kalman filter
• Bayesian estimation via adaptive RWMH (bayes_dsge()) with optional parallel chain execution (n_cores)
• Particle filter and PMMH (bayes_particle()) for fully nonlinear Bayesian estimation without linearization
• Bayes factor model comparison (bayes_factor()) with Kass-Raftery evidence scales and posterior model probabilities
• Variance decomposition (variance_decomposition()) -- both unconditional steady-state shares and forecast-error variance decomposition at user-supplied horizons
• Optimal Simple Rules (osr()) -- finds the optimal coefficients in a user-specified, restricted policy rule (e.g. Taylor rule coefficients) given a quadratic welfare loss
• Conditional forecasts (conditional_forecast()) -- forecasts conditional on a pre-specified path for a subset of observables (Waggoner-Zha 1999 minimum-norm shocks)
• IRF matching estimation (irf_match()) -- estimate structural parameters by matching DSGE impulse responses to a user-supplied target (e.g. VAR-estimated IRFs), Christiano-Eichenbaum-Evans style
• DSGE-VAR (bayes_dsge_var(), bayes_dsge_var_mh()) -- Bayesian VAR with DSGE-implied prior (Del Negro & Schorfheide 2004). The MH variant jointly estimates structural parameters, shock SDs, and the prior weight lambda; forecast() and conditional_forecast() methods support fan-chart projections and path-conditioned forecasts -- Dynare-parity workflow
• Anticipated / news shocks in perfect_foresight_nonlinear()
• Derived parameters in dsge_model(derived = ...) for encoding models with Dynare-style # macro substitutions (e.g. Smets-Wouters 2007), with the doubling-algorithm Lyapunov solver handling highly-persistent state dynamics robustly
• Calibrated smoother (calibrated_smoother()) -- Kalman smoothing on a model that has been solved but not estimated
• predetermined() -- alias for state-variable declarations, matching Dynare's predetermined_variables
• Perfect foresight with expectation errors (perfect_foresight_expect_err()) -- realised vs. subjective shock paths with revising expectations
• Extended path simulation (extended_path()) -- stochastic simulation by solving a perfect-foresight path per period
• Endogenous priors (endogenous_prior()) -- Christiano-Trabandt- Walentin prior on model-implied second moments
• Global sensitivity analysis (global_sensitivity()) -- Sobol' indices and Morris elementary effects
• Discretionary optimal policy (discretionary_policy()) -- time-consistent (no-commitment) feedback rule via the Soederlind / Dennis fixed-point iteration
• Generalised/Simulated Method of Moments (gmm_estimate(), smm_estimate()) -- moment-matching estimators with optional two-step weighting
• Sequential Monte Carlo sampler (bayes_smc()) -- tempered SMC for robust posterior sampling, especially for multimodal posteriors
• LaTeX model export (model_latex()) -- write a model's equations to compilable LaTeX, with Greek substitution and time subscripts
• Skewed Kalman filter (kalman_filter_skewed()) -- exact third-cumulant propagation for skew-normal structural shocks
• Markov-switching volatility (ms_filter()) -- regime-switching shock variances via the Kim (1994) filter, with smoothed regime probabilities
• PAC equations (pac_weights(), pac_target_loading(), pac_simulate()) -- FRB/US-style polynomial adjustment costs with closed-form forward-looking weights
• Ramsey optimal policy via linear-quadratic regulator (ramsey_policy(), welfare_loss())
• Second- and third-order perturbation with pruned simulation
• Occasionally binding constraints (OccBin) for ZLB and other bounds
• Perfect foresight deterministic transition paths -- both linearized (perfect_foresight()) and fully nonlinear via stacked-time Newton (perfect_foresight_nonlinear())
• Kalman smoothing and historical shock decomposition
• Local identification diagnostics and parameter sensitivity analysis
• Robust (sandwich) standard errors for ML estimation
• Posterior predictive checks and marginal likelihood
• Model-implied covariance matrices and prediction tools
• Publication-ready plots with a unified theme: forecast fan charts with history, +/- 2 sigma bands on smoothed states, and overlays for comparing perfect-foresight paths

Installation

# Install from CRAN
install.packages("dsge")

# Or install the development version from GitHub
# install.packages("devtools")
devtools::install_github("Mustapha-Wasseja/dsge-package")

Quick Start: Maximum Likelihood

library(dsge)

# Define a simple New Keynesian model
nk <- dsge_model(
  obs(p   ~ beta * lead(p) + kappa * x),       # Phillips curve
  unobs(x ~ lead(x) - (r - lead(p) - g)),      # IS curve
  obs(r   ~ psi * p + u),                       # Taylor rule
  state(u ~ rhou * u),                          # Monetary shock
  state(g ~ rhog * g),                          # Demand shock
  fixed = list(beta = 0.99),
  start = list(kappa = 0.1, psi = 1.5, rhou = 0.7, rhog = 0.9)
)

# Estimate by maximum likelihood
fit <- estimate(nk, data = your_data)
summary(fit)

# Postestimation
irf(fit, periods = 20) |> plot()               # Impulse responses
forecast(fit, horizon = 12) |> plot()           # Forecasts
smooth_states(fit) |> plot()                    # Kalman-smoothed states
shock_decomposition(fit) |> plot()              # Historical decomposition
check_identification(fit)                       # Identification diagnostics
robust_vcov(fit)                                # Sandwich standard errors
model_covariance(fit)                           # Model-implied moments

Bayesian Estimation

# Specify priors
my_priors <- list(
  kappa = prior("beta", shape1 = 30, shape2 = 70),
  psi   = prior("gamma", shape = 184, rate = 122.7),
  rhou  = prior("beta", shape1 = 70, shape2 = 20),
  rhog  = prior("beta", shape1 = 70, shape2 = 20)
)

# Run MCMC
fit_bayes <- bayes_dsge(nk, data = your_data, priors = my_priors,
                        chains = 2, iter = 10000, warmup = 5000)

# Diagnostics and results
summary(fit_bayes)
plot(fit_bayes, type = "trace")
plot(fit_bayes, type = "prior_posterior")
plot(fit_bayes, type = "irf", periods = 20)

Supported priors: normal, beta, gamma, uniform, inv_gamma.

Nonlinear DSGE Models

rbc <- dsgenl_model(
  "1/C = beta / C(+1) * (alpha * exp(Z) * K^(alpha-1) + 1 - delta)",
  "K(+1) = exp(Z) * K^alpha - C + (1 - delta) * K",
  "Z(+1) = rho * Z",
  observed = "C",
  endo_state = "K",
  exo_state = "Z",
  fixed = list(alpha = 0.33, beta = 0.99, delta = 0.025),
  start = list(rho = 0.9),
  ss_guess = c(C = 2, K = 30, Z = 0)
)

sol <- solve_dsge(rbc, params = c(alpha = 0.33, beta = 0.99,
                                   delta = 0.025, rho = 0.9),
                  shock_sd = c(Z = 0.01))
irf(sol, periods = 40) |> plot()

Advanced Features

Second- and Third-Order Perturbation

sol2 <- solve_dsge(rbc, params = params, shock_sd = sd, order = 2)
simulate_2nd_order(sol2, periods = 200)
irf_2nd_order(sol2, periods = 40)

sol3 <- solve_dsge(rbc, params = params, shock_sd = sd, order = 3)
simulate_3rd_order(sol3, periods = 200)

Occasionally Binding Constraints

# ZLB constraint on the interest rate
obc <- simulate_occbin(sol,
  constraints = list("r >= 0"),
  shocks = list(g = -0.05),
  horizon = 40)
plot(obc)

Perfect Foresight Paths

# Linearized perfect foresight (fast, small shocks)
pf <- perfect_foresight(sol,
  shocks = list(Z = c(-0.05, -0.03, -0.01)),
  horizon = 60)
plot(pf)

# Fully nonlinear perfect foresight via stacked-time Newton
# (recommended for large shocks where nonlinearities matter)
pf_nl <- perfect_foresight_nonlinear(rbc,
  params   = c(rho = 0.9),
  shock_sd = c(Z = 0.01),
  shocks   = list(Z = 0.10),
  horizon  = 40)
plot(pf_nl)

Particle Filter and PMMH

# Bootstrap particle filter likelihood for nonlinear models
ll <- particle_filter_loglik(sol2, data = your_data, n_particles = 1000)

# Particle Marginal Metropolis-Hastings -- fully nonlinear Bayesian
fit_pmmh <- bayes_particle(rbc, data = your_data, priors = my_priors,
                           n_particles = 500, chains = 2, iter = 5000)

Ramsey Optimal Policy

# Quadratic welfare loss on inflation and output gap
rp <- ramsey_policy(sol,
  Q_xx = diag(c(p = 1, x = 0.5)),
  Q_yy = diag(c(r = 0.1)))
welfare_loss(sol, rp$F)   # evaluate welfare under the optimal rule

Bayes Factor Model Comparison

# Compare two Bayesian fits
bf <- bayes_factor(fit_bayes_A, fit_bayes_B,
                   prior_odds = c(0.5, 0.5))
print(bf)   # log Bayes factor + Kass-Raftery evidence label

Variance Decomposition

# Unconditional decomposition (long-run shares)
vd <- variance_decomposition(sol)
print(vd)
plot(vd)

# Forecast-error variance decomposition at multiple horizons
fevd <- variance_decomposition(sol, horizon = c(1, 4, 8, 20))
plot(fevd)

Optimal Simple Rules

# Optimal Taylor-rule coefficient
res <- osr(nk,
  params      = c(kappa = 0.1, psi = 1.5, rhou = 0.7, rhog = 0.9),
  shock_sd    = c(e.u = 1.0, e.g = 0.5),
  osr_params  = c(psi = 1.5),
  welfare_weights = list(Q_yy = c(p = 1, x = 0.5, r = 0.1)),
  lower = 1.01, upper = 5.0)
print(res)

Conditional Forecasts

# Hold the policy rate at 0.5 for the next 4 periods
cf <- conditional_forecast(fit, horizon = 12,
  condition = list(r = c(0.5, 0.5, 0.5, 0.5, rep(NA, 8))))
plot(cf)

IRF Matching Estimation

# Match the DSGE IRF to an externally estimated target
est <- irf_match(nk,
  params_start   = c(kappa = 0.15, psi = 2.0, rhou = 0.6, rhog = 0.8),
  shock_sd_start = c(e.u = 0.8, e.g = 0.7),
  target         = target_irf_dataframe)

DSGE-VAR

# (a) Conditional-on-(theta,lambda) Bayesian VAR with DSGE prior
fit_dv <- bayes_dsge_var(sol, data = your_data,
                         p = 4, lambda = 1.0, n_draws = 1000)
print(fit_dv)
# Compare different lambdas by marginal likelihood
sapply(c(0.5, 1, 2, 5),
       function(l) bayes_dsge_var(sol, your_data, p=4, lambda=l,
                                  n_draws=200)$log_marg_lik)

# (b) Joint MH estimation of (theta_DSGE, sigma, lambda) -- Dynare-parity
priors <- list(kappa = prior("beta",  shape1 = 2, shape2 = 8),
               psi   = prior("normal", mean = 1.5, sd = 0.25),
               rhou  = prior("beta",  shape1 = 2, shape2 = 2),
               rhog  = prior("beta",  shape1 = 8, shape2 = 2))
fit_mh <- bayes_dsge_var_mh(nk, data = your_data, priors = priors,
                            p = 4, chains = 2, iter = 2000)
print(fit_mh)

# Unconditional + conditional forecasts from the DSGE-VAR posterior
fc_unc  <- forecast(fit_mh, horizon = 12)
fc_cond <- conditional_forecast(fit_mh, horizon = 12,
            condition = list(r = c(0.5, 0.5, 0.5, 0.5, rep(NA, 8))))
plot(fc_unc); plot(fc_cond)

Anticipated / News Shocks (Perfect Foresight)

# Nonlinear PF correctly anticipates a future TFP shock
pf_news <- perfect_foresight_nonlinear(rbc,
  params   = c(rho = 0.9),
  shock_sd = c(Z = 0.01),
  shocks   = list(Z = c(0, 0, 0.05)),  # announced at t=1, arrives at t=3
  horizon  = 40)
plot(pf_news)

Importing Dynare Models

# Read a Dynare .mod file: model, calibration, steady state, shocks, priors
rbc <- read_dynare(system.file("examples", "rbc.mod", package = "dsge"))
rbc                        # summary, auxiliary variables and any notes

sol <- solve_dsge(rbc)     # solves at the file's calibration
plot(irf(sol, periods = 40))

# Files with estimated_params / varobs estimate directly, using the
# translated priors (data columns named as in Dynare)
# fit <- bayes_dsge(read_dynare("model.mod"), data = my_data)

# Optimal policy and occasionally binding constraints declared in the file
# ram <- solve_dsge(read_dynare("ramsey.mod"))       # ramsey_model
# opt <- osr(read_dynare("osr.mod"))                  # osr_params, optim_weights
# zlb <- simulate_occbin(read_dynare("zlb.mod"))      # occbin_constraints

Dynare timing is handled automatically (lags become auxiliary state variables, so k(-1) capital timing needs no rewriting), and macro directives (@#define, @#for, @#if, ...) are expanded in R. MATLAB statements in the file (calibrations, verbatim blocks) and MATLAB steady-state files (<model>_steadystate.m, with fsolve and helper functions) are run by a built-in MATLAB interpreter. Results have been checked against Dynare 6.0 (see dev/dynare-validation/): the first-order impulse responses of every model in Johannes Pfeifer's DSGE_mod collection that Dynare runs in Octave, second- and third-order decision rules of 21 nonlinear models, the Smets-Wouters (2007) likelihood (with lik_init = 1 and 2), and a full Bayesian estimation.

Parallel MCMC Chains

# Run chains in parallel across cores (PSOCK on Windows, fork on POSIX)
fit_bayes <- bayes_dsge(nk, data = your_data, priors = my_priors,
                        chains = 4, iter = 10000, warmup = 5000,
                        n_cores = 4)

Feature Comparison

Documentation

• CRAN: https://CRAN.R-project.org/package=dsge
• Vignette: vignette("introduction", package = "dsge")
• Examples: system.file("examples", package = "dsge")

References

• Andrieu, C., Doucet, A. and Holenstein, R. (2010). "Particle Markov chain Monte Carlo methods." Journal of the Royal Statistical Society: Series B, 72(3), 269-342.
• Gordon, N. J., Salmond, D. J. and Smith, A. F. M. (1993). "Novel approach to nonlinear/non-Gaussian Bayesian state estimation." IEE Proceedings F, 140(2), 107-113.
• Juillard, M., Laxton, D., McAdam, P. and Pioro, H. (1998). "An algorithm competition: First-order iterations versus Newton-based techniques." Journal of Economic Dynamics and Control, 22, 1291-1318.
• Kass, R. E. and Raftery, A. E. (1995). "Bayes factors." Journal of the American Statistical Association, 90(430), 773-795.
• Klein, P. (2000). "Using the generalized Schur form to solve a multivariate linear rational expectations model." Journal of Economic Dynamics and Control, 24(10), 1405-1423.
• Schmitt-Grohe, S. and Uribe, M. (2004). "Solving dynamic general equilibrium models using a second-order approximation to the policy function." Journal of Economic Dynamics and Control, 28(4), 755-775.

License

MIT