## Reproduce p1-fit through p1-sim-manual outside knit to surface the error.
suppressPackageStartupMessages({
  library(tidyverse); library(lubridate); library(MASS); library(lmtest)
})

approval <- read_csv("data/approval.csv", show_col_types = FALSE)
cat("approval colnames:\n"); print(colnames(approval))

df <- approval |>
  mutate(date        = ym(paste(year, month, sep = "-")),
         oil         = avg.price,
         approve_lag = lag(approve, 1),
         oil_lag     = lag(oil, 1)) |>
  drop_na()
cat("df rows:", nrow(df), "\n")

fit_ardl <- lm(approve ~ approve_lag + oil + oil_lag +
                 sept.oct.2001 + iraq.war,
               data = df)
cat("fit_ardl coefnames:\n"); print(names(coef(fit_ardl)))

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))
cat("df_cf rows:", nrow(df_cf), "\n")

anchor_idx <- which(df_cf$date < T_star)
fc_idx     <- which(df_cf$date >= T_star)
H          <- length(fc_idx)
cat("anchor n:", length(anchor_idx), " fc n:", H, "\n")

set.seed(2025); S <- 50
beta_hat <- coef(fit_ardl); V_hat <- vcov(fit_ardl)
draws    <- mvrnorm(n = S, mu = beta_hat, Sigma = V_hat)
cat("draws colnames:\n"); print(colnames(draws))

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)

cat("attempt s=1, i=fc_idx[1]\n")
s <- 1; i <- fc_idx[1]; b <- draws[s, ]
val <-
  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("val =", val, "  class:", class(val), "\n")
