inst/doc/qpmR-estimation.R

In qpmR: Quarterly Projection Models for Monetary Policy Analysis

## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>",
                      fig.width = 7, fig.height = 5)
set.seed(1)

## -----------------------------------------------------------------------------
library(qpmR)
pr <- priors(
  rho = beta(0.5, 0.2),
  e   = invgamma(1, 0.5)    # a shock name means that shock's sd
)
pr

## -----------------------------------------------------------------------------
m0 <- qpm_model(variables = vars(x = "x"), shocks = shocks(e),
                equations = eqs(x ~ rho * x[-1] + e),
                params = list(rho = 0.5))
m_true <- qpm_calibrate(m0, rho = 0.8, sigma = c(e = 1.5))
obs <- simulate(qpm_solve(m_true), nsim = 250, seed = 4)

est <- qpm_estimate(m0, obs, pr, iter = 800, chains = 2, seed = 5,
                    verbose = FALSE)
est

## -----------------------------------------------------------------------------
plot(est)

## -----------------------------------------------------------------------------
m_hat <- apply_estimate(est, "mean")
round(coef(est, "mean"), 3)

## -----------------------------------------------------------------------------
fc <- posterior_forecast(est, horizon = 10, ndraws = 80)
plot(fc, vars = "x")

## -----------------------------------------------------------------------------
m_bad <- qpm_model(variables = vars(x = "x"), shocks = shocks(e),
                   equations = eqs(x ~ a * b * x[-1] + e),
                   params = list(a = 0.6, b = 0.9))
qpm_identify(m_bad, params = c("a", "b"))

## -----------------------------------------------------------------------------
qpm_identify(qpm_template("bkl"),
             params = c("b1", "b2", "b3", "c1", "c2", "a1", "a3"),
             observables = c("pi", "i", "q", "y_gap", "dy_obs"))

## -----------------------------------------------------------------------------
ml_ar1 <- marginal_likelihood(est)
ml_ar1

# a deliberately misspecified rival: white noise (rho fixed at 0)
est_wn <- qpm_estimate(qpm_calibrate(m0, rho = 0), obs,
                       priors(e = invgamma(1, 0.5)),
                       iter = 800, chains = 2, seed = 6, verbose = FALSE)
ml_wn <- marginal_likelihood(est_wn)
cat(sprintf("log Bayes factor, AR(1) vs white noise: %.1f\n",
            ml_ar1$logml - ml_wn$logml))

## ----eval = FALSE-------------------------------------------------------------
# mcz <- qpm_calibrate(qpm_template("bkl", trends = "rw"),
#                      pi_tar = 2, istar_ss = 2, pistar_ss = 2,
#                      prem_ss = 1, a5 = 0.4)
# cz <- czechia[czechia$period >= "1999",
#               c("period", "pi4", "i", "q", "dy_obs", "istar", "pistar")]
# est_cz <- 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, 0.5)
# ), iter = 3000, chains = 2, seed = 42)