R/dynare.R

In qpmR: Quarterly Projection Models for Monetary Policy Analysis

#' Export a model to a Dynare .mod file
#'
#' Writes the model as Dynare source so that anyone can re-solve it in
#' Dynare and compare. This is deliberate: qpmR's solver is pinned in its
#' own tests to analytic solutions, but "does it agree with Dynare?" is
#' the first question a modeller asks, and the answer should be cheap to
#' check rather than something to take on trust.
#'
#' The original equations are exported, not qpmR's internal first-order
#' system: Dynare creates its own auxiliary variables for long leads and
#' lags, so the two implementations agree only if both handle them
#' correctly. Expectation wrappers are dropped (`E(pi[+1])` becomes
#' `pi(+1)`), since leads in Dynare are already model-consistent
#' expectations.
#'
#' @param model A `qpm_model`.
#' @param file Output path. The default `NULL` returns the Dynare source as
#'   a character vector; a file is written only when a path is given.
#' @param irf Horizon for `stoch_simul`; `0` omits impulse responses.
#' @param order Approximation order passed to `stoch_simul` (the models
#'   are linear, so first order is exact).
#' @param extra Character vector of extra Dynare statements appended
#'   after `stoch_simul`, e.g. code to export results.
#' @return The file path invisibly, or the Dynare source as a character
#'   vector when `file` is `NULL`.
#' @examples
#' src <- write_dynare(qpm_template("bkl"), file = NULL)
#' cat(head(src, 15), sep = "\n")
#' @export
write_dynare <- function(model, file = NULL, irf = 40, order = 1,
                         extra = NULL) {
  stopifnot(inherits(model, "qpm_model"))
  nms <- c(model$vars$name, model$shocks, names(model$params))
  bad <- nms[!grepl("^[A-Za-z_][A-Za-z0-9_]*$", nms)]
  if (length(bad))
    stop(sprintf(paste0("these names are not valid Dynare identifiers: %s\n",
                        "  Dynare allows letters, digits and underscores only (no dots)."),
                 paste(bad, collapse = ", ")), call. = FALSE)

  wrapv <- function(x, width = 72) {
    out <- character(0); cur <- ""
    for (w in x) {
      cand <- if (nzchar(cur)) paste(cur, w) else w
      if (nchar(cand) > width) { out <- c(out, cur); cur <- w } else cur <- cand
    }
    c(out, cur)
  }

  lines <- c(
    "// Generated by qpmR::write_dynare()",
    sprintf("// model: %s", model$name),
    sprintf("// qpmR %s", as.character(utils::packageVersion("qpmR"))),
    "//",
    "// Cross-check: solve this in Dynare and compare its IRFs with",
    "// qpmR::irf() on the same model. They should agree to solver tolerance.",
    "",
    "var", wrapv(model$vars$name), ";", "",
    "varexo", wrapv(model$shocks), ";", ""
  )

  if (length(model$params)) {
    lines <- c(lines, "parameters", wrapv(names(model$params)), ";", "")
    lines <- c(lines,
               sprintf("%s = %s;", names(model$params),
                       vapply(unname(model$params), fmt_g, character(1))), "")
  }

  lines <- c(lines, "model(linear);")
  for (f in model$equations)
    lines <- c(lines, paste0("  ", dynare_eq_text(f)))
  lines <- c(lines, "end;", "")

  lines <- c(lines, "shocks;",
             sprintf("  var %s; stderr %s;", model$shocks,
                     vapply(unname(model$sigma[model$shocks]), fmt_g, character(1))),
             "end;", "")

  lines <- c(lines, "steady;", "check;", "")
  so <- sprintf("stoch_simul(order=%d%s, nograph, noprint);", order,
                if (irf > 0) sprintf(", irf=%d", irf) else ", irf=0")
  lines <- c(lines, so)
  if (!is.null(extra)) lines <- c(lines, "", extra)

  if (is.null(file)) return(lines)
  dir.create(dirname(normalizePath(file, mustWork = FALSE)),
             recursive = TRUE, showWarnings = FALSE)
  writeLines(lines, file)
  invisible(file)
}

fmt_g <- function(x) trimws(format(x, digits = 15, scientific = FALSE))

# Rewrite one qpmR equation as Dynare source: x[-1] -> x(-1), x[+1] -> x(+1),
# E() dropped. The AST is transformed and then deparsed, so R's own operator
# precedence rules produce correctly parenthesised output.
dynare_eq_text <- function(f) {
  rec <- function(e) {
    if (!is.call(e)) return(e)
    op <- e[[1L]]
    opname <- if (is.symbol(op)) as.character(op) else ""
    if (opname == "[") {
      nm <- as.character(e[[2L]])
      # numeric, not integer: an integer literal would deparse as "-1L"
      k <- as.numeric(eval(e[[3L]], envir = baseenv()))
      if (k == 0) return(as.symbol(nm))
      return(as.call(list(as.symbol(nm), k)))
    }
    if (opname == "E") return(rec(e[[2L]]))
    out <- e
    for (i in seq_along(e)[-1L]) out[[i]] <- rec(e[[i]])
    out
  }
  txt <- paste0(deparse1(rec(f[[2L]])), " = ", deparse1(rec(f[[3L]])), ";")
  # deparse writes a one-quarter lead as name(1); make it explicit for Dynare
  gsub("([A-Za-z_][A-Za-z0-9_]*)\\(([0-9]+)\\)", "\\1(+\\2)", txt)
}