Skip to contents

qpmR

Quarterly Projection Models for Monetary Policy Analysis in R.

qpmR builds, solves, filters, estimates and reports the semi-structural quarterly projection models (QPM) used in central-bank Forecasting and Policy Analysis Systems (FPAS): output gap, Phillips curve, forward-looking policy rule, and exchange-rate block, solved under model-consistent expectations.

The goal is the workflow, not just the equations:

data -> filtering -> gaps -> model -> baseline -> judgment -> scenarios -> report

Every stage of that chain is implemented, exercised end to end on a real quarterly dataset for Czechia that ships with the package, and cross-checked against Dynare. Documentation: https://mustapha-wasseja.github.io/qpmR/.

Installation 

install.packages("qpmR")            # from CRAN, once accepted

# development version
# install.packages("pak")
pak::pak("Mustapha-Wasseja/qpmR")

Quickstart 

library(qpmR)

m <- qpm_template("bkl")          # canonical small open economy QPM
summary(m)

sol <- qpm_solve(m)               # QZ solution + Blanchard-Kahn check
sol

ir <- irf(sol, shock = "eps_i")   # 100bp-style policy tightening
plot(ir, vars = c("pi", "y_gap", "i", "q"))

histq <- simulate(sol, nsim = 48, seed = 7, burn = 20)
fc <- qpm_forecast(sol, from = histq, horizon = 12)
plot(fc, vars = c("pi", "i", "y_gap", "q"))

# transmission experiment: double exchange-rate pass-through
m2 <- qpm_calibrate(m, b3 = 0.2)
plot(irf(qpm_solve(m2), shock = "eps_q"), vars = c("pi", "i"))

The vignettes walk through the whole chain: vignette("qpmR-quickstart") takes a forecast round from data to report, and vignette("qpmR-estimation") covers priors, posteriors, identification and Bayes factors.

What it does

The model layer

  • A native R model language. Declare gap-form models with qpm_model(): lags/leads as x[-1] / E(x[+1]), labels and units on every variable, calibration with typo protection (qpm_calibrate()). Long lags and leads (e.g. E(pi4[+4]) in the policy rule) are handled automatically via auxiliary states.
  • A generalized-Schur (QZ) solver (Klein 2000) with honest Blanchard–Kahn diagnostics: indeterminacy and explosiveness are reported as typed errors that say which economics usually causes them (Taylor principle, pure-UIP unit roots). Unit roots are counted explicitly, so random-walk trends solve.
  • The canonical Berg–Karam–Laxton small-open-economy QPM as a calibrated template: qpm_template("bkl") — IS curve, hybrid Phillips curve, forward-looking inflation-targeting rule, dampened UIP, and trend/foreign processes, for an illustrative emerging-economy calibration (“Meridia”, 5% inflation target). trends = "rw" makes the equilibrium exchange rate and potential growth unit-root processes.
  • Model dynamics: irf() with peak-effect tables, fevd() for forecast error variance decompositions, model_properties() for model-implied standard deviations and autocorrelations set against the data’s, simulate(), qpm_forecast() with analytic fan bands, steady_state(), eigen_table(), and qpm_lint() for specification checks.

Filtering

  • qpm_filter() — Kalman filter + RTS smoother over the solved model: jointly infers every latent state (output gap, neutral rate, equilibrium exchange rate, trends) and the historical structural shocks from whatever subset of variables you observe. Missing data and ragged edges are handled; innovation diagnostics (Ljung–Box, outlier flags) are printed. Unit-root models switch to a diffuse prior automatically. The likelihood is tested against the exact closed-form Gaussian likelihood.
  • qpm_decompose() — exact historical shock decompositions of the smoothed history (contributions verified to sum to the states), with stacked-bar charts.
  • state_space() — the exact T, R, Z, H, Qc, P1 matrices used internally, exported so other estimators can build on qpmR.
  • residuals(), fitted(), logLik() and nobs() on a filtration, so AIC() and BIC() work without special handling.
  • A real country dataset — czechia, quarterly from 1996 in model units, compiled reproducibly from FRED/OECD/Eurostat/ECB public endpoints. Filtering it reproduces the known history: the pre-GFC boom, the 2009 and 2013 recessions, the COVID crater, the koruna’s trend real appreciation, and the post-GFC fall in potential growth — all pinned in the test suite.
mcz <- qpm_calibrate(qpm_template("bkl", trends = "rw"),
                     pi_tar = 2, istar_ss = 2, pistar_ss = 2, prem_ss = 1)
cz <- czechia[czechia$period >= "1999",
              c("period", "pi4", "i", "q", "dy_obs", "istar", "pistar")]
fit <- qpm_filter(mcz, cz)
plot(fit, vars = c("y_gap", "dy_bar", "r_bar", "q_gap"))
plot(qpm_decompose(fit), var = "pi4")

Forecasting and policy analysis

  • qpm_condition() — hard conditional forecasts with the minimum-norm implied shocks reported in standard deviations, an explicit anticipated switch (announced-at-start vs period-by-period surprises — materially different in any forward-looking model), instrument restrictions, and conditional fan bands that collapse at conditioned points. The anticipation recursion is verified against brute-force perfect foresight in the test suite.
  • qpm_scenario() — shock-based alternatives, announced or surprise.
  • add_judgment() / judgment_log() — judgment as a first-class, logged operation: state the change, qpmR back-solves the supporting shocks, records author/time/rationale, and flags anything requiring more than two standard deviations.
  • qpm_risk() / risk_log() — a balance of risks: bands become two-piece normal, so the mode stays on the projection while the mean shifts by the stated skew with total variance held fixed.
  • Forecast rounds — qpm_round() archives the whole pipeline (model + calibration + data vintage + filtration + conditioned forecast) as one replayable object; save_round() / load_round() / list_rounds() manage a plain-directory store with CSV sidecars for auditing without R.
  • compare_rounds() — the revision decomposition. “Inflation for 2027-Q1 is 0.55pp higher than we said in March: +0.15 new outturns, +0.40 judgment.” Computed by re-running the pipeline swapping one ingredient at a time (parameters → data revisions → new data → conditions → judgment); the contributions telescope exactly and the endpoints are verified against the archived rounds.
base <- qpm_forecast(qpm_solve(mcz), from = fit, horizon = 12)
hold <- qpm_condition(base, i = c("2026-Q3" = 3.5, "2026-Q4" = 3.5),
                      anticipated = TRUE, instruments = "eps_i")

rA <- qpm_round("2026-Q1 March", mcz, cz_to_2025Q4, horizon = 12)
rB <- add_judgment(qpm_round("2026-Q3 September", mcz, cz, horizon = 12),
                   pi4 = c("2027-Q1" = 0.4), author = "prices desk",
                   rationale = "announced energy-tariff increase")
compare_rounds(rA, rB)

Estimation

  • priors() — a prior mini-language (beta, gamma, invgamma, normal, uniform, truncate) in the mean/sd parametrization economists write down, scoped inside priors() so base R’s beta() and gamma() are never masked.
  • qpm_estimate() — Bayesian estimation of any subset of parameters and shock sds over the Kalman-filter likelihood: posterior mode, adaptive random-walk Metropolis seeded by the BFGS Hessian, split R-hat / ESS diagnostics, and a “learned” column comparing posterior to prior spread. method = "mle" uses the same machinery. coef(), vcov(), confint(), summary() and logLik() work on the result.
  • posterior_forecast() — fan charts that integrate over the posterior: every draw re-solves the model and re-filters the data.
  • qpm_identify() — Iskrev-style identification diagnostics before sampling: which parameters have no effect, which are only jointly identified, where the Jacobian loses rank.
  • marginal_likelihood() — modified harmonic mean with a Laplace cross-check; differences across models are log Bayes factors.
est <- qpm_estimate(mcz, cz, priors(
  b1 = beta(0.70, 0.10), b2 = gamma(0.25, 0.10), b3 = gamma(0.10, 0.05),
  c1 = beta(0.70, 0.10), c2 = truncate(normal(1.5, 0.25), lower = 1),
  eps_pi = invgamma(1, 0.5)
), iter = 4000, chains = 2)
est
plot(est)                                   # prior vs posterior
plot(posterior_forecast(est, horizon = 12)) # parameter-uncertainty fans

Country adaptation

  • add_block() — extension blocks that adapt a template to a country without forking it, composable and recorded in the model.
  • block_food_cpi() — headline CPI split into food and core, the configuration every LIC/EM engagement needs and rebuilds by hand. A food supply shock moves headline while leaving core untouched.
  • block_fx_intervention() — a managed float: one model spanning free float through peg by a single intensity argument.
  • qpm_diff() — a country customization reviewed as a diff of structure; qpm_compare_models() compares behaviour (impulse responses and implied moments).
  • qpm_disaggregate() — Denton-Cholette and Chow-Lin temporal disaggregation of annual data to quarterly, for the many economies that publish national accounts only annually.
m <- add_block(qpm_template("bkl"), block_food_cpi(weight = 0.45))
plot(irf(qpm_solve(m), shock = "eps_pifood"),
     vars = c("pi_food", "pi", "pi_core", "i"))
qpm_diff(qpm_template("bkl"), m)

Reporting, audit and policy experiments

  • qpm_report() — a round becomes the monetary policy report: an executive summary with the numbers filled in, fan charts, gaps, the shock decomposition, the judgment ledger, the revision against the previous round, and a reproducibility appendix. The .Rmd source is always written so teams edit the text, not the plumbing; rendering degrades gracefully where pandoc is unavailable.
  • chart_pack() — the standard round chart set as a multi-page PDF.
  • verify_round() — re-runs an archived round and confirms the published numbers still come back, flagging version drift and hand-edited CSV sidecars.
  • qpm_rule_eval() — score alternative policy rules over a grid by the unconditional loss, computed exactly from the stationary covariance, and trace the inflation-output variability frontier.
  • qpm_counterfactual() — replay history with shocks switched off or scaled: “what if the central bank had simply followed its rule?”
  • write_dynare() — export any model as a Dynare .mod file.
verify_round(r)                       # does the archive still reproduce?
chart_pack(r, "chart_pack.pdf")
qpm_report(r, "mpr.html", compare_to = previous_round)

Does it agree with Dynare?

Yes — and you can check rather than take it on trust. write_dynare() exports any model as a .mod file, and the test suite pins qpmR’s impulse responses to golden files produced by Dynare 6.0:

write_dynare(qpm_template("bkl"), "bkl.mod")

Across 4080 impulse-response points (12 shocks x 17 variables x 20 quarters) the largest discrepancy is 1.4e-14, with steady states agreeing to 8e-14. The food-block model agrees to 1.7e-14, which also validates add_block() independently — Dynare parses the original equations and builds its own auxiliary variables, so the two implementations agree only if both handle long leads and lags correctly. Regenerate the golden files with data-raw/dynare_golden.R.

Speed

The Kalman filter and the stationary-covariance solve are compiled (RcppArmadillo). Both keep reference implementations in R, and the test suite pins the compiled paths to them to machine precision — that agreement is what makes the fast path trustworthy. Use options(qpmR.use_cpp = FALSE) to fall back to R.

A posterior draw on the Czech model (22 states, 110 quarters) costs 27 ms, against 181 ms before this work — so a 6000-draw estimate takes under three minutes rather than eighteen. Two of the three gains were not the filter:

before after
model assembly (build_first_order) 27.5 ms 1.7 ms
stationary covariance (Lyapunov) 39 ms 1 ms
Kalman filter 55 ms 16 ms

The structure of the first-order system does not depend on the parameters, so it is computed once per model and cached; and the Lyapunov equation is solved by O(N^3) squaring rather than the O(N^6) Kronecker system, which also speeds up model_properties(), qpm_rule_eval() and qpm_identify().

Testing

More than 800 assertions across 28 files, at about 90% line coverage — measured on every push by the test-coverage workflow. The suite pins the solver to analytic solutions (AR(1), hybrid roots, brute-force perfect foresight), the Kalman filter to the exact closed-form Gaussian likelihood and the compiled filter to its R reference, the revision decomposition to exact telescoping, and the whole solver to Dynare (above).

Design commitments

  1. Everything has an escape hatch. sol$P, sol$Q, eigen_table() and state_space() expose the actual matrices; nothing is hidden in closures.
  2. Errors teach. A Blanchard–Kahn failure names the economics that usually causes it, not just the rank condition it violates.
  3. Verification is cheap. The test suite pins the solver to analytic solutions and to Dynare for the shipped templates, and every compiled path to its R reference implementation.

Roadmap

Version Focus
0.1 Model DSL, QZ solver, BK diagnostics, IRFs, simulation, forecasts, BKL template
0.2 Kalman filter/smoother, shock decompositions, unit-root trends with diffuse initialization, real country dataset (czechia)
0.3 Conditional forecasts (anticipated vs unanticipated), scenarios, judgment ledger, forecast rounds, round store, revision decomposition
0.4 Bayesian estimation (priors, adaptive RWM, R-hat/ESS), identification diagnostics, marginal likelihood, posterior fans, estimation vignette
1.0 Country adaptation (extension blocks, qpm_diff()), reporting (qpm_report(), chart_pack()) and audit (verify_round())
1.1 Compiled Kalman filter and Lyapunov solver, fevd(), model_properties(), balance of risks, rule evaluation, counterfactuals, model comparison, temporal disaggregation, standard R generics, Dynare cross-check, pkgdown site; CRAN submission
next Exact Durbin-Koopman diffuse initialization

References

License

MIT.