## Full p1-sim-manual under same lib-load order as the Rmd setup chunk.
suppressPackageStartupMessages({
  library(tidyverse); library(lubridate)
  library(plm); library(fixest); library(nlme); library(lme4)
  library(parameters); library(sandwich); library(lmtest); library(tseries)
  library(MASS); library(simcf)
  library(ggdag); library(dagitty)
  library(patchwork); library(cowplot); library(kableExtra); library(texreg)
})

approval <- read_csv("data/approval.csv", show_col_types = FALSE)
df <- approval |>
  mutate(date        = ym(paste(year, month, sep = "-")),
         oil         = avg.price,
         approve_lag = lag(approve, 1),
         oil_lag     = lag(oil, 1)) |>
  drop_na()
fit_ardl <- lm(approve ~ approve_lag + oil + oil_lag +
                 sept.oct.2001 + iraq.war, data = df)
T_star    <- ymd("2003-03-01")
oil_base  <- mean(df$oil[df$date < T_star])
df_cf <- df |>
  mutate(oil_cf     = if_else(date >= T_star, oil_base, oil),
         oil_cf_lag = lag(oil_cf, 1)) |>
  mutate(oil_cf_lag = if_else(is.na(oil_cf_lag), oil_lag, oil_cf_lag))

anchor_idx <- which(df_cf$date < T_star)
fc_idx     <- which(df_cf$date >= T_star)
H          <- length(fc_idx)

# Predicted-value step
approve_cf <- df_cf$approve
beta <- coef(fit_ardl)
for (i in fc_idx) {
  newrow <- tibble(
    approve_lag   = approve_cf[i - 1],
    oil           = df_cf$oil_cf[i],
    oil_lag       = df_cf$oil_cf_lag[i],
    sept.oct.2001 = df_cf$sept.oct.2001[i],
    iraq.war      = df_cf$iraq.war[i]
  )
  approve_cf[i] <- predict(fit_ardl, newdata = newrow)
}
cat("predict-loop OK, last forecast =", tail(approve_cf, 1), "\n")

# Manual MASS sim, full S=1000
set.seed(2025); S <- 1000
beta_hat <- coef(fit_ardl); V_hat <- vcov(fit_ardl)
draws    <- mvrnorm(n = S, mu = beta_hat, Sigma = V_hat)

sim_path   <- matrix(NA_real_, nrow = S, ncol = nrow(df_cf))
sim_path[, anchor_idx] <- matrix(df_cf$approve[anchor_idx],
                                 nrow = S, ncol = length(anchor_idx),
                                 byrow = TRUE)
for (s in seq_len(S)) {
  b <- draws[s, ]
  for (i in fc_idx) {
    sim_path[s, i] <-
      b["(Intercept)"]       +
      b["approve_lag"]   * sim_path[s, i - 1] +
      b["oil"]           * df_cf$oil_cf[i]    +
      b["oil_lag"]       * df_cf$oil_cf_lag[i] +
      b["sept.oct.2001"] * df_cf$sept.oct.2001[i] +
      b["iraq.war"]      * df_cf$iraq.war[i]
  }
}
cat("sim_path ok; mean of last col =", mean(sim_path[, ncol(sim_path)]), "\n")

manual_cf <- tibble(
  date     = df_cf$date,
  observed = df_cf$approve,
  pe       = apply(sim_path, 2, mean),
  lo       = apply(sim_path, 2, quantile, probs = 0.025),
  hi       = apply(sim_path, 2, quantile, probs = 0.975)
)
print(tail(manual_cf, 3))
