## Reproduce p1-sim-simcf to get the actual error message.
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)
})
lag <- dplyr::lag; filter <- dplyr::filter; select <- dplyr::select

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)

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)

cat("draws colnames:\n"); print(colnames(draws))
cat("Inspect cfMake signature:\n"); print(args(cfMake))
cat("Inspect ldvsimev signature:\n"); print(args(ldvsimev))

formula_cf <- approve ~ oil + oil_lag + sept.oct.2001 + iraq.war
xhyp <- cfMake(formula_cf, df_cf, H)
cat("xhyp class:", class(xhyp), " names:\n"); print(names(xhyp))
cat("xhyp$x dim:", dim(xhyp$x), "\n")
cat("xhyp$x colnames:\n"); print(colnames(xhyp$x))

for (h in seq_len(H)) {
  xhyp <- cfChange(xhyp, "oil",           x = df_cf$oil_cf[fc_idx[h]],     scen = h)
  xhyp <- cfChange(xhyp, "oil_lag",       x = df_cf$oil_cf_lag[fc_idx[h]], scen = h)
  xhyp <- cfChange(xhyp, "sept.oct.2001", x = df_cf$sept.oct.2001[fc_idx[h]], scen = h)
  xhyp <- cfChange(xhyp, "iraq.war",      x = df_cf$iraq.war[fc_idx[h]],   scen = h)
}
cat("xhyp$x[1,]:\n"); print(xhyp$x[1, ])

beta_for_sim <- draws[, c("(Intercept)", "oil", "oil_lag",
                          "sept.oct.2001", "iraq.war")]
phi_hat   <- coef(fit_ardl)["approve_lag"]
y_anchor  <- df_cf$approve[max(anchor_idx)]
cat("phi_hat =", phi_hat, " y_anchor =", y_anchor, "\n")

sim_simcf <- ldvsimev(
  xhyp     = xhyp,
  simbetas = beta_for_sim,
  ci       = 0.95,
  constant = 1,
  phi      = phi_hat,
  lagY     = y_anchor
)
cat("sim_simcf names:\n"); print(names(sim_simcf))
cat("pe[1:5]:\n"); print(sim_simcf$pe[1:5])
