## Reproduce the exact failing sequence with prints around the mvrnorm call.
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)

cat("BEFORE predict-loop:\n")
cat("  df.residual:", df.residual(fit_ardl), "\n")
cat("  sigma:      ", sigma(fit_ardl), "\n")
cat("  vcov diag:  ", diag(vcov(fit_ardl)), "\n\n")

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)

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("AFTER predict-loop:\n")
cat("  df.residual:", df.residual(fit_ardl), "\n")
cat("  sigma:      ", sigma(fit_ardl), "\n")
cat("  vcov diag:  ", diag(vcov(fit_ardl)), "\n\n")

cat("identical(beta, coef(fit_ardl))? ", identical(beta, coef(fit_ardl)), "\n")
cat("Are beta values printed:\n"); print(beta)
cat("Was 'beta' overwritten anywhere? class(beta) =", class(beta), "\n")
