rm(list = ls()); set.seed(123)

# =============================
# B1: Baseline VAR -> NGML SVAR
# =============================

suppressPackageStartupMessages({ library(vars); library(svars); library(zoo) })

# --- Paths (WD = project_root/Code) ---
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 data ---
D  <- read.csv(file.path(in_dir, "kaenzig_baseline_data.csv"),
               header = TRUE, stringsAsFactors = FALSE, na.strings = c("NaN","NA",""))
ym <- as.yearmon(gsub("M","-", D$time), "%Y-%m")

Y  <- D[, 2:7]
Y[] <- lapply(Y, as.numeric)
ok <- complete.cases(Y)
Y  <- Y[ok, , drop = FALSE]
ym <- ym[ok]

Y_ts <- ts(as.matrix(Y),
           start = c(as.integer(format(ym[1], "%Y")), as.integer(format(ym[1], "%m"))),
           frequency = 12)

proxy_full   <- scan(file.path(in_dir, "kaenzig_baseline_proxy.csv"), quiet = TRUE)
news_full    <- read.csv(file.path(in_dir, "kaenzig_baseline_structural_shock.csv"),
                         header = TRUE, stringsAsFactors = FALSE)$shock
timeNum_full <- read.csv(file.path(in_dir, "kaenzig_baseline_hd.csv"), header = FALSE)[[1]]

# --- VAR + NGML SVAR ---
p_lags  <- 12
n_ahead <- 50
nboot   <- 10000

plain.var <- VAR(Y_ts, p = p_lags, type = "const")
u <- resid(plain.var)
T <- nrow(u)

proxy   <- tail(proxy_full,   T)
news    <- tail(news_full,    T)
timeNum <- tail(timeNum_full, T)

ngml <- id.ngml(plain.var)
B <- ngml$B
K <- ncol(B)

eps <- t(solve(B, t(u)))
colnames(eps) <- paste0("eps_", seq_len(ncol(eps)))

# --- df + se ---
wcsv(data.frame(shock = colnames(eps),
                df    = as.numeric(ngml$df),
                df_se = as.numeric(ngml$df_SE)),
     "thispaper_baseline_ngml_df_and_se.csv")

# --- Bootstrap ---
boot <- mb.boot(ngml, design = "recursive", b.length = 24,
                n.ahead = n_ahead, nboot = nboot, nc = 1,
                dd = NULL, signrest = NULL,
                itermax = 300, steptol = 200, iter2 = 50)

# --- Oil news shock (max |corr| with Kaenzig news) ---
r_news    <- apply(eps, 2, function(x) cor(x, news))
oil_j     <- which.max(abs(r_news))
oil_shock <- as.numeric(eps[, oil_j])

# --- IRFs (own responses) + Känzig-style median-recentered percentile bands ---
true_df <- as.data.frame(boot$true$irf)
h1 <- suppressWarnings(as.numeric(true_df[[1]]))
has_h <- (ncol(true_df) >= 2) && all(!is.na(h1)) && h1[1] == 0
h <- if (has_h) h1 else 0:(nrow(true_df) - 1)

irf_true <- as.matrix(if (has_h) true_df[, -1, drop = FALSE] else true_df)
own_idx  <- seq(1, ncol(irf_true), by = K)
IRF0     <- irf_true[, own_idx, drop = FALSE]

H  <- nrow(IRF0)
kk <- ncol(IRF0)

nb  <- length(boot$bootstrap)
IRFb <- array(NA_real_, c(H, kk, nb))

for (b in seq_len(nb)) {
  draw <- as.data.frame(boot$bootstrap[[b]]$irf)
  
  hd <- suppressWarnings(as.numeric(draw[[1]]))
  has_h_draw <- (ncol(draw) >= 2) && all(!is.na(hd)) && hd[1] == 0
  if (has_h_draw) draw <- draw[, -1, drop = FALSE]
  
  M <- as.matrix(draw)[, own_idx, drop = FALSE]
  if (nrow(M) != H) {
    if (nrow(M) > H) M <- M[1:H, , drop = FALSE]
    else M <- rbind(M, matrix(NA_real_, H - nrow(M), ncol(M)))
  }
  IRFb[, , b] <- M
}

# --- quantiles recentered by bootstrap median ---
Q <- function(p) {
  out <- matrix(NA_real_, H, kk)
  for (i in seq_len(H)) for (j in seq_len(kk))
    out[i, j] <- as.numeric(quantile(IRFb[i, j, ], probs = p, na.rm = TRUE))
  out
}

MED <- Q(0.50)

# 90% (alpha = 0.10): [0.05, 0.95]
LO90_raw <- Q(0.05)
HI90_raw <- Q(0.95)
LO90 <- LO90_raw - MED + IRF0
HI90 <- HI90_raw - MED + IRF0

# 68% (alpha = 0.32): [0.16, 0.84]
LO68_raw <- Q(0.16)
HI68_raw <- Q(0.84)
LO68 <- LO68_raw - MED + IRF0
HI68 <- HI68_raw - MED + IRF0

wcsv(data.frame(h = h),   "thispaper_baseline_irf_horizon.csv")
wcsv(as.data.frame(IRF0), "thispaper_baseline_irf_point.csv")
wcsv(as.data.frame(LO90), "thispaper_baseline_irf_lo90.csv")
wcsv(as.data.frame(HI90), "thispaper_baseline_irf_hi90.csv")
wcsv(as.data.frame(LO68), "thispaper_baseline_irf_lo68.csv")
wcsv(as.data.frame(HI68), "thispaper_baseline_irf_hi68.csv")

# --- Historical decomposition (point) ---
hd0 <- hd(ngml, series = 1, transition = 0)$hidec
wcsv(data.frame(date            = as.Date(as.yearmon(time(hd0))),
                actual_demeaned = as.numeric(hd0[, 1]),
                hd_point        = as.numeric(hd0[, 2])),
     "thispaper_baseline_hd_point.csv")

# --- Core exports ---
wcsv(data.frame(timeNum = timeNum, structural_shock = oil_shock),
     "thispaper_baseline_structural_shock.csv")

wcsv(cbind(timeNum = timeNum, as.data.frame(eps)),
     "thispaper_baseline_structural_shocks.csv")

wcsv(as.data.frame(B), "thispaper_baseline_B.csv")

wcsv(data.frame(shock = colnames(eps),
                corr_with_kaenzig_news = as.numeric(r_news)),
     "thispaper_baseline_corr_with_kaenzig_news.csv")

wcsv(data.frame(shock = colnames(eps),
                corr_with_proxy = apply(eps, 2, function(x) cor(x, proxy))),
     "thispaper_baseline_corr_with_proxy.csv")

wcsv(as.data.frame(u), "thispaper_baseline_VAR_residuals.csv")

cat("Done. CSV outputs in: ", out_dir, "\n", sep = "")