# NGML SVAR (baseline 6 + 1 extension) and max-|corr| with Kaenzig extension shock

rm(list = ls()); set.seed(123)
suppressPackageStartupMessages({ library(vars); library(svars); library(zoo) })

# --- Paths ---
project_dir <- dirname(getwd())
in_dir  <- file.path(project_dir, "Data", "kaenzig_exports")
out_dir <- file.path(project_dir, "Data", "thispaper_exports")
dir.create(out_dir, showWarnings = FALSE, recursive = TRUE)

wcsv <- function(x, fn) write.csv(x, file.path(out_dir, fn), row.names = FALSE)

# --- Load + merge data ---
Db <- read.csv(file.path(in_dir, "kaenzig_baseline_data.csv"),
               header = TRUE, stringsAsFactors = FALSE, na.strings = c("NaN","NA",""))
De <- read.csv(file.path(in_dir, "kaenzig_extensions_data.csv"),
               header = TRUE, stringsAsFactors = FALSE, na.strings = c("NaN","NA",""))

base_vars <- names(Db)[2:7]
ext_vars  <- setdiff(names(De), "time")

M <- merge(Db[, c("time", base_vars)], De, by = "time", all = FALSE, sort = FALSE)
ym_all <- as.yearmon(gsub("M","-", M$time), "%Y-%m")

# --- Settings ---
p_lags <- 12
rows <- vector("list", length(ext_vars))

# --- Loop over extensions ---
for (i in seq_along(ext_vars)) {
  
  x <- ext_vars[i]
  
  Ydf <- M[, c(base_vars, x)]
  Ydf[] <- lapply(Ydf, as.numeric)
  
  first_ok <- which(!is.na(Ydf[[x]]))[1]
  Ydf <- Ydf[first_ok:nrow(Ydf), , drop = FALSE]
  ym  <- ym_all[first_ok:length(ym_all)]
  
  ok  <- complete.cases(Ydf)
  Ydf <- Ydf[ok, , drop = FALSE]
  ym  <- ym[ok]
  
  Y_ts <- ts(as.matrix(Ydf),
             start = c(as.integer(format(ym[1], "%Y")),
                       as.integer(format(ym[1], "%m"))),
             frequency = 12)
  
  plain.var <- VAR(Y_ts, p = p_lags, type = "const")
  u <- resid(plain.var)
  T <- nrow(u)
  
  B <- id.ngml(plain.var)$B
  eps <- t(solve(B, t(u)))
  
  ka_full <- read.csv(
    file.path(in_dir, paste0("kaenzig_extensions_", x, "_structural_shock.csv"))
  )$shock
  
  ka <- tail(as.numeric(ka_full), T)
  
  r <- apply(eps, 2, function(s) cor(s, ka, use = "pairwise.complete.obs"))
  rows[[i]] <- data.frame(
    extension = x,
    corr_with_kaenzig = abs(r[which.max(abs(r))])
  )
}

# --- Export ---
wcsv(do.call(rbind, rows), "thispaper_extensions_shock_correlations.csv")