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

## ----simulate-----------------------------------------------------------------
set.seed(42)
true_rho <- 0.8
true_sd <- 1.0
n <- 200

z <- numeric(n)
for (i in 2:n) z[i] <- true_rho * z[i - 1] + true_sd * rnorm(1)
dat <- data.frame(y = z)

## ----specify------------------------------------------------------------------
library(dsge)

m <- dsge_model(
  obs(y ~ z),
  state(z ~ rho * z),
  start = list(rho = 0.5)
)

print(m)

## ----solve--------------------------------------------------------------------
sol <- solve_dsge(m, params = c(rho = 0.8))
print(sol)

## ----estimate-----------------------------------------------------------------
fit <- estimate(m, data = dat)
summary(fit)

## ----postest------------------------------------------------------------------
# Policy matrix (maps states to controls)
policy_matrix(fit, se = FALSE)

# Transition matrix (state dynamics)
transition_matrix(fit, se = FALSE)

# Stability check
stability(fit)

## ----irf, fig.height=3--------------------------------------------------------
irfs <- irf(fit, periods = 20, se = FALSE)
plot(irfs)

## ----forecast, fig.height=3---------------------------------------------------
fc <- forecast(fit, horizon = 12)
plot(fc)

## ----nk-model-----------------------------------------------------------------
nk <- dsge_model(
  obs(p   ~ beta * lead(p) + kappa * x),
  unobs(x ~ lead(x) - (r - lead(p) - g)),
  obs(r   ~ psi * p + u),
  state(u ~ rhou * u),
  state(g ~ rhog * g),
  fixed = list(beta = 0.96),
  start = list(kappa = 0.1, psi = 1.5, rhou = 0.7, rhog = 0.9)
)

print(nk)

## ----nk-sim-est---------------------------------------------------------------
# Solve at true parameters to get state-space form
true_params <- c(beta = 0.96, kappa = 0.085, psi = 1.94,
                 rhou = 0.70, rhog = 0.95)
true_sd <- c(u = 2.3, g = 0.57)
sol <- solve_dsge(nk, params = true_params, shock_sd = true_sd)

# Simulate data
set.seed(123)
n <- 200
states <- matrix(0, n, 2)
colnames(states) <- c("u", "g")
for (t in 2:n) {
  states[t, "u"] <- 0.70 * states[t - 1, "u"] + 2.3 * rnorm(1)
  states[t, "g"] <- 0.95 * states[t - 1, "g"] + 0.57 * rnorm(1)
}
Z <- sol$D %*% sol$G
obs_data <- states %*% t(Z)
colnames(obs_data) <- c("p", "r")
nk_dat <- as.data.frame(obs_data)

# Estimate
nk_fit <- estimate(nk, data = nk_dat,
                   start = list(kappa = 0.1, psi = 1.5,
                                rhou = 0.7, rhog = 0.9),
                   control = list(maxit = 500))
summary(nk_fit)

## ----nk-irf, fig.height=5-----------------------------------------------------
nk_irfs <- irf(nk_fit, periods = 20)
plot(nk_irfs)

## ----nk-forecast, fig.height=4------------------------------------------------
nk_fc <- forecast(nk_fit, horizon = 12)
plot(nk_fc)

## ----rbc-model----------------------------------------------------------------
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)
)
print(rbc)

## ----rbc-ss-------------------------------------------------------------------
params <- c(alpha = 0.33, beta = 0.99, delta = 0.025, rho = 0.9)
ss <- steady_state(rbc, params = params)
print(ss)

## ----rbc-solve----------------------------------------------------------------
sol <- solve_dsge(rbc, params = params, shock_sd = c(Z = 0.01))
print(sol)
stability(sol)

## ----rbc-irf, fig.height=3----------------------------------------------------
rbc_irfs <- irf(sol, periods = 40, se = FALSE)
plot(rbc_irfs)

## ----bayes-priors, eval=FALSE-------------------------------------------------
# my_priors <- list(
#   beta  = prior("beta", shape1 = 95, shape2 = 5),
#   kappa = prior("beta", shape1 = 30, shape2 = 70),
#   psi   = prior("gamma", shape = 20, rate = 13.3),
#   rhou  = prior("beta", shape1 = 70, shape2 = 20),
#   rhog  = prior("beta", shape1 = 70, shape2 = 20)
#   # shock SDs default to inv_gamma(0.01, 0.01)
# )

## ----bayes-estimate, eval=FALSE-----------------------------------------------
# fit_bayes <- bayes_dsge(nk, data = nk_dat, priors = my_priors,
#                         chains = 2, iter = 10000, warmup = 5000,
#                         seed = 42)
# 
# summary(fit_bayes)  # posterior table with ESS, R-hat, MCSE

## ----bayes-diag, eval=FALSE---------------------------------------------------
# # MCMC diagnostics
# plot(fit_bayes, type = "trace")            # trace plots (all chains)
# plot(fit_bayes, type = "density")          # posterior density + prior overlay
# plot(fit_bayes, type = "prior_posterior")   # dedicated prior-vs-posterior comparison
# plot(fit_bayes, type = "running_mean")     # cumulative mean convergence
# plot(fit_bayes, type = "acf")              # autocorrelation diagnostics
# plot(fit_bayes, type = "pairs")            # pairwise scatter + correlations
# plot(fit_bayes, type = "all")              # combined trace + density panel
# 
# # Parameter selection (works with all types above)
# plot(fit_bayes, type = "trace", pars = c("kappa", "psi"))
# plot(fit_bayes, type = "density", pars = "rhou")
# 
# # Posterior IRFs with credible bands
# plot(fit_bayes, type = "irf", periods = 12, n_draws = 500)
# # Or equivalently:
# bayes_irfs <- irf(fit_bayes, periods = 12, n_draws = 500)
# plot(bayes_irfs)