# ============================================================ # # I_LP_baseline.R # # Script I — LP Baseline # # Project: Macroeconomic vs Financial Uncertainty, Euro Area # # # # MAIN SPECIFICATION: # # MU shock = ar_meu_innov (Comunale-Nguyen AR innov) # # FU uncertainty = fu_uncertainty_innov (BH GARCH-based) # # FU risk premium = fu_riskpremium_innov (BH VRP) # # # # ROBUSTNESS (PCA_CES): # # MU shock = MEU_PCA_CES_innov (PCA: MEU+EPU+CESIEUR) # # FU shocks unchanged # # # # ============================================================ # # ── Auto working directory (RStudio only) ───────────────────── # When shared with supervisor: opens script → Session auto-set if (requireNamespace("rstudioapi", quietly = TRUE) && rstudioapi::isAvailable()) { script_path <- rstudioapi::getSourceEditorContext()$path if (nchar(script_path) > 0) setwd(dirname(script_path)) } # ── Packages ────────────────────────────────────────────────── required_pkgs <- c("tidyverse", "readxl", "openxlsx", "lubridate", "sandwich", "lmtest", "patchwork", "AER") new_pkgs <- required_pkgs[!required_pkgs %in% installed.packages()[,"Package"]] if (length(new_pkgs)) install.packages(new_pkgs) suppressPackageStartupMessages({ library(tidyverse) library(readxl) library(openxlsx) library(lubridate) library(sandwich) library(lmtest) library(patchwork) library(AER) }) # ── Namespace conflict resolution ───────────────────────────── select <- dplyr::select filter <- dplyr::filter summarise <- dplyr::summarise lag <- dplyr::lag recode <- dplyr::recode # ── Output folders ──────────────────────────────────────────── FIG_MAIN <- "output/figures/A1_baseline" FIG_PCA <- "output/figures/PCA_CESI_robustness" TAB_MAIN <- "output/tables/A1_baseline" TAB_PCA <- "output/tables/PCA_CESI_robustness" for (d in c(FIG_MAIN, FIG_PCA, TAB_MAIN, TAB_PCA)) dir.create(d, recursive = TRUE, showWarnings = FALSE) # Helper: print AND save ggplot save_fig <- function(plot, filename, dir = FIG_MAIN, width = 10, height = 7) { print(plot) ggsave(file.path(dir, filename), plot, width = width, height = height) cat("Saved:", file.path(dir, filename), "\n") } # ── LaTeX output folder ─────────────────────────────────────── LAT_MAIN <- "output/latex/A1_baseline" dir.create(LAT_MAIN, recursive = TRUE, showWarnings = FALSE) # Helper: save xlsx AND a LaTeX-ready CSV alongside it # Usage: save_table(df, "my_table") # Produces: output/tables/A1_baseline/my_table.xlsx # output/latex/A1_baseline/my_table.tex (kable booktabs) save_table <- function(df, name, dir = TAB_MAIN, latex_dir = LAT_MAIN, caption = NULL, label = NULL, digits = 4) { # 1. Excel as before openxlsx::write.xlsx(df, file.path(dir, paste0(name, ".xlsx")), overwrite = TRUE) # 2. LaTeX via knitr::kable if (!requireNamespace("knitr", quietly = TRUE)) install.packages("knitr") if (!requireNamespace("kableExtra", quietly = TRUE)) install.packages("kableExtra") cap <- if (!is.null(caption)) caption else name lab <- if (!is.null(label)) label else paste0("tab:", gsub("[^a-zA-Z0-9]", "_", name)) tex <- knitr::kable( df, format = "latex", booktabs = TRUE, digits = digits, caption = cap, label = lab ) %>% kableExtra::kable_styling(latex_options = c("hold_position", "scale_down")) writeLines(as.character(tex), file.path(latex_dir, paste0(name, ".tex"))) cat("Saved:", file.path(dir, paste0(name, ".xlsx")), "\n") cat("Saved:", file.path(latex_dir, paste0(name, ".tex")), "\n") } # ============================================================ # # 0. LOAD ANALYSIS DATASET # # Single source of truth built by 00_build_dataset.R # # ============================================================ # cat("=======================================================\n") cat(" I_LP_BASELINE — Euro Area Uncertainty LP\n") cat("=======================================================\n\n") cat("--- 0. Loading analysis dataset ---\n") dat <- read_csv("data/final/analysis_dataset.csv", show_col_types = FALSE) %>% mutate(date = as.Date(date)) %>% arrange(date) cat(sprintf("Dataset: %d obs x %d vars | %s to %s\n\n", nrow(dat), ncol(dat), format(min(dat$date), "%Y-%m"), format(max(dat$date), "%Y-%m"))) # ============================================================ # # 1. ANALYSIS SAMPLE # # 2003M6–2024M7: start determined by MEU availability # # Drop remaining NAs in all core LP variables # # ============================================================ # cat("--- 1. Trimming to analysis sample ---\n") lp_data <- dat %>% filter(date >= as.Date("2003-06-01"), date <= as.Date("2024-07-01")) %>% arrange(date) cat(sprintf("LP sample before outcome-specific NA handling: %d obs | %s to %s\n\n", nrow(lp_data), format(min(lp_data$date), "%Y-%m"), format(max(lp_data$date), "%Y-%m"))) # Quick completeness check on all variables used below cat("Missing values in LP variables:\n") lp_vars_check <- c("ar_meu_innov", "fu_uncertainty_innov", "fu_riskpremium_innov", "MEU_PCA_CES_innov", "ip_yoy", "d_inflation", "d_unemployment", "capital_goods_growth", "nfc_growth", "d_policy_rate", "output_gap_std", "ciss_std", "d_credit_spreads", "covid_dummy") lp_data %>% select(any_of(lp_vars_check)) %>% summarise(across(everything(), ~sum(is.na(.)))) %>% pivot_longer(everything(), names_to = "variable", values_to = "n_missing") %>% filter(n_missing > 0) %>% print() cat("(If empty: no missing values in key LP variables)\n\n") # ============================================================ # # UNCERTAINTY MEASURES CORRELATION DIAGNOSTIC # # ============================================================ # cat("=== Correlation between uncertainty measures ===\n\n") innov_vars <- c( "ar_meu_innov", # MU: Communale-Nguyen AR shock "MEU_PCA_CES_innov", # MU: PCA-CES AR shock "fu_uncertainty_innov", # FU level: BH expected physical variance "fu_riskpremium_innov" # FU vol: BH variance risk premium (C-S) ) # Keep only variables that exist in dataset innov_vars <- intersect(innov_vars, names(lp_data)) cor_mat <- lp_data %>% select(all_of(innov_vars)) %>% drop_na() %>% cor() cat("Pairwise correlations:\n") print(round(cor_mat, 3)) # Flag any concerning correlations cat("\nCorrelation flags:\n") upper_tri <- cor_mat[upper.tri(cor_mat)] pairs <- which(upper.tri(cor_mat), arr.ind = TRUE) for (i in seq_len(nrow(pairs))) { r <- cor_mat[pairs[i, 1], pairs[i, 2]] v1 <- rownames(cor_mat)[pairs[i, 1]] v2 <- colnames(cor_mat)[pairs[i, 2]] if (abs(r) > 0.40) { cat(sprintf(" |r| > 0.40 — %s ~ %s: %.3f (potential overlap)\n", v1, v2, r)) } } cat(" (|r| < 0.20 = well-separated | 0.20-0.40 = acceptable | > 0.40 = review)\n\n") # Save write.xlsx( as.data.frame(round(cor_mat, 3)), "output/tables/A1_baseline/shock_correlations.xlsx", overwrite = TRUE ) # ============================================================ # # 2. SHOCK DEFINITIONS # # ============================================================ # # All shocks are standardized innovation series. # They remove predictable own dynamics from the corresponding # uncertainty measure, but they are NOT necessarily mutually # orthogonal. Therefore, all three shocks enter jointly in the LP. # # MAIN SPECIFICATION: # ar_meu_innov — Comunale-Nguyen MEU shock # AR(2) w/ macro controls, LB p=0.126 # fu_uncertainty_innov — BH expected physical variance innovation. # (FU level) This captures unexpected movements in expected financial # market variance under the physical probability measure. AR(2), LB p=0.126 # fu_riskpremium_innov — BH variance risk premium innovation. # (FU vol, Caballero-Simsek channel) # This captures unexpected movements in the excess of # risk-neutral implied variance over expected physical variance # AR(6), LB p=0.155 # ==================================================================== # SHOCKS_MAIN <- c("ar_meu_innov", "fu_uncertainty_innov", "fu_riskpremium_innov") SHOCKS_PCA <- c("MEU_PCA_CES_innov", "fu_uncertainty_innov", "fu_riskpremium_innov") # Display names for plots shock_labels <- c( ar_meu_innov = "MU shock (MEU innov.)", MEU_PCA_CES_innov = "MU shock (PCA-CES innov.)", fu_uncertainty_innov = "FU uncertainty (BH expected var.)", fu_riskpremium_innov = "FU risk premium (BH VRP)" ) shock_colors <- c( ar_meu_innov = "#2166ac", MEU_PCA_CES_innov = "#4393c3", fu_uncertainty_innov = "#d6604d", fu_riskpremium_innov = "#4dac26" ) # ============================================================ # # 3. CONTROL SPECIFICATION # # # # Economic rationale for each control: # # # # d_policy_rate (lagged): # # Monetary policy controls for the ECB reaction function. # # First difference used: policy rate is I(1) in sample. # # Bernanke & Blinder (1992): rate changes propagate to # # investment and credit. Enter lagged: policy works with # # a lag and is predetermined relative to current shock. # # # # output_gap_std (lagged): # # Business cycle position. Uncertainty shocks may have # # state-dependent effects (Bloom 2009, Caballero-Simsek). # # Standardised: allows comparison across specifications. # # Enter lagged: output gap is pre-determined at t. # # # # ciss_std (lagged): # # Composite Indicator of Systemic Stress (Hollo et al 2012).# # Controls for broad financial conditions separately from # # the BH FU measures. CISS enters lagged — its current # # value may be endogenous to the FU shocks. # # # # d_credit_spreads (lagged): # # Credit premium (Gilchrist & Zakrajšek 2012 spirit). # # Captures tightening in credit conditions beyond CISS. # # First difference: credit_spreads appears non-stationary. # # # # nflation (lagged): # # Price dynamics. Controls for supply-driven inflation that # # may co-move with uncertainty (energy price shocks). # # Not included as control for the inflation LP itself # # # # covid_dummy (contemporaneous): # # Exogenous structural break: March–May 2020. # # The three COVID months are outliers in every series; # # including a dummy prevents them from dominating the IRF. # # Enters contemporaneously: the shock IS the dummy event. # # ============================================================ # # Baseline state controls — these enter with lags only BASE_CONTROLS <- c( "d_policy_rate", "output_gap_std", "ciss_std", "d_credit_spreads", "inflation" ) # Optional sentiment control, if available ESI_CONTROL <- intersect(c("esi_std", "d_esi", "esi"), names(lp_data)) # Contemporaneous exogenous dummy only EXOG <- c("covid_dummy") # LP settings H <- 24 # N_LAGS will be selected below using VAR lag-selection diagnostics. # If you want to force a manual lag length, set MANUAL_N_LAGS to an integer. MANUAL_N_LAGS <- NA_integer_ MAX_LAG_SELECT <- 12 # Baseline uses 6 monthly lags as a conservative dynamic specification. # Earlier lag-selection exercises suggested shorter lag lengths, so # lag-length sensitivity should be reported separately. # Outcome-specific control function get_controls <- function(outcome, spec = c("baseline", "extended")) { spec <- match.arg(spec) controls <- BASE_CONTROLS # Add lagged real activity where economically required if (outcome %in% c("inflation", "unemployment", "capital_goods_growth", "nfc_growth", "ciss_std", "ciss", "credit_spreads", "credit_spreads_std", "d_credit_spreads")) { controls <- c(controls, "ip_yoy") } # If outcome is inflation, remove inflation from controls if (outcome %in% c("inflation", "d_inflation")) { controls <- setdiff(controls, c("inflation", "d_inflation")) } # If outcome is CISS, remove CISS from controls if (outcome %in% c("ciss", "ciss_std", "d_ciss")) { controls <- setdiff(controls, c("ciss", "ciss_std", "d_ciss")) } # If outcome is credit spreads, remove credit-spread controls if (outcome %in% c("credit_spreads", "credit_spreads_std", "d_credit_spreads")) { controls <- setdiff( controls, c("credit_spreads", "credit_spreads_std", "d_credit_spreads") ) } # If outcome is IP, no need to add ip_yoy as control because # run_lp already adds own lags of the outcome. if (outcome == "ip_yoy") { controls <- setdiff(controls, "ip_yoy") } # Extended sentiment specification if (spec == "extended" && length(ESI_CONTROL) > 0) { controls <- c(controls, ESI_CONTROL[1]) } # Keep only variables that actually exist controls <- intersect(unique(controls), names(lp_data)) return(controls) } cat("--- 3. LP specification ---\n") cat("Shocks (main): ", paste(SHOCKS_MAIN, collapse=", "), "\n") cat("Shocks (pca-ces): ", paste(SHOCKS_PCA, collapse=", "), "\n") cat("Base controls: ", paste(BASE_CONTROLS, collapse=", "), "\n") cat("ESI control: ", ifelse(length(ESI_CONTROL) > 0, ESI_CONTROL[1], "not available"), "\n") cat("Contemporaneous: ", paste(EXOG, collapse=", "), "\n") cat("Horizon H: ", H, "\n") cat("Lag length: ", ifelse(is.na(MANUAL_N_LAGS), "to be selected (VAR below)", as.character(MANUAL_N_LAGS)), "\n\n") # ============================================================ # # 3B. REQUIRED VARIABLE CHECKS # # ============================================================ # # Run this before the OUTCOMES block to diagnose cat("CISS candidates in data:", intersect(c("ciss_std","ciss","d_ciss"), names(lp_data)), "\n") cat("Credit candidates in data:", intersect(c("credit_spreads_std","credit_spreads","d_credit_spreads"), names(lp_data)), "\n") cat("Other outcome candidates:", intersect(c("ip_yoy","inflation","unemployment", "capital_goods_growth","nfc_growth"), names(lp_data)), "\n") # Decide available CISS and credit-spread outcomes CISS_OUTCOME <- if ("ciss_std" %in% names(lp_data)) { "ciss_std" } else if ("ciss" %in% names(lp_data)) { "ciss" } else { NA_character_ } CREDIT_OUTCOME <- if ("credit_spreads_std" %in% names(lp_data)) { "credit_spreads_std" } else if ("credit_spreads" %in% names(lp_data)) { "credit_spreads" } else if ("d_credit_spreads" %in% names(lp_data)) { "d_credit_spreads" } else { NA_character_ } OUTCOMES <- c( "ip_yoy", "inflation", "unemployment", "capital_goods_growth", "nfc_growth", CISS_OUTCOME, CREDIT_OUTCOME ) %>% purrr::discard(is.na) %>% intersect(names(lp_data)) required_vars <- unique(c( SHOCKS_MAIN, SHOCKS_PCA, OUTCOMES, BASE_CONTROLS, EXOG )) missing_required <- setdiff(required_vars, names(lp_data)) if (length(missing_required) > 0) { stop("Missing required variables: ", paste(missing_required, collapse = ", ")) } cat("Required-variable check passed.\n\n") cat("Missing values in required LP variables:\n") lp_data %>% select(all_of(required_vars)) %>% summarise(across(everything(), ~ sum(is.na(.)))) %>% pivot_longer(everything(), names_to = "variable", values_to = "n_missing") %>% filter(n_missing > 0) %>% print(n = Inf) cat("(If empty: no missing values in required variables)\n\n") # ============================================================ # # 3C. LP LAG-LENGTH SELECTION # # ============================================================ # # Purpose: # Choose the number of lags used inside run_lp(): # - lags of the outcome # - lags of controls # # We use a compact monthly macro-financial state vector to avoid # overfitting. The selection is diagnostic; robustness checks should # compare shorter and longer lag choices. # ============================================================ # cat("\n=== 3C. LP LAG-LENGTH SELECTION ===\n") lag_select_vars <- c( "ip_yoy", "inflation", "unemployment", "d_policy_rate", "output_gap_std", "ciss_std", "d_credit_spreads" ) lag_select_vars <- intersect(lag_select_vars, names(lp_data)) lag_select_df <- lp_data %>% select(all_of(lag_select_vars)) %>% drop_na() cat("Lag-selection variables:\n") print(lag_select_vars) cat(sprintf("Lag-selection sample: %d observations\n", nrow(lag_select_df))) if (nrow(lag_select_df) < 80) { warning("Lag-selection sample is small. Interpret VAR lag selection carefully.") } lag_sel <- vars::VARselect( lag_select_df, lag.max = MAX_LAG_SELECT, type = "const" ) cat("\nVAR lag-selection recommendations:\n") print(lag_sel$selection) cat("\nFull lag-selection criteria:\n") print(round(lag_sel$criteria, 3)) N_LAGS_AIC <- as.integer(lag_sel$selection["AIC(n)"]) N_LAGS_HQ <- as.integer(lag_sel$selection["HQ(n)"]) N_LAGS_BIC <- as.integer(lag_sel$selection["SC(n)"]) N_LAGS_FPE <- as.integer(lag_sel$selection["FPE(n)"]) lag_decision_table <- tibble( criterion = c("AIC", "HQ", "BIC/SC", "FPE"), selected_lag = c(N_LAGS_AIC, N_LAGS_HQ, N_LAGS_BIC, N_LAGS_FPE) ) cat("\nCompact lag decision table:\n") print(lag_decision_table) # ------------------------------------------------------------ # # Baseline rule # ------------------------------------------------------------ # # BIC/SC is usually preferred for parsimony in finite samples. # AIC often selects more lags. For monthly macro LPs, 3–6 lags # is usually a reasonable range, but we do not want to impose 6 # without checking. # # Recommended baseline: # - Use BIC if BIC >= 2. # - If BIC selects 1, use max(2, HQ) to avoid an overly thin # dynamic structure. # - Always report robustness with 6 lags. # ------------------------------------------------------------ # if (!is.na(MANUAL_N_LAGS)) { N_LAGS <- MANUAL_N_LAGS lag_choice_reason <- paste0("manual choice: ", MANUAL_N_LAGS) } else { N_LAGS <- max(2, N_LAGS_BIC) lag_choice_reason <- paste0( "baseline rule: max(2, BIC-selected lag = ", N_LAGS_BIC, ")" ) } # Robustness lag set LAG_ROBUSTNESS_SET <- sort(unique(c( N_LAGS, N_LAGS_HQ, N_LAGS_AIC, 6 ))) # Avoid absurdly large AIC choices LAG_ROBUSTNESS_SET <- LAG_ROBUSTNESS_SET[LAG_ROBUSTNESS_SET <= 8] cat("\nChosen baseline N_LAGS:", N_LAGS, "\n") cat("Reason:", lag_choice_reason, "\n") cat("Lag robustness set:", paste(LAG_ROBUSTNESS_SET, collapse = ", "), "\n\n") write.xlsx( list( selection = lag_decision_table, criteria = as.data.frame(round(lag_sel$criteria, 3)) ), file.path(TAB_MAIN, "lp_lag_selection.xlsx"), overwrite = TRUE ) # ============================================================ # # 4. INSTRUMENT FIRST-STAGE DIAGNOSTICS # # Tests whether available instruments are sufficiently # # correlated with each shock for LP-IV identification. # # # This is a screening diagnostic only. The formal LP-IV script # should replicate the exact LP-IV first stage horizon by horizon. # # H0: instrument coefficient = 0 (instrument irrelevant) # # H1: instrument coefficient ≠ 0 (instrument relevant) # # # # Strength thresholds (Staiger-Stock 1997): # # F > 16.38 : strong (10% Stock-Yogo size) # # F > 10 : acceptable rule of thumb # # F < 10 : weak — interpret IV results with caution # # # # We report both OLS F and HAC-robust F. # # With persistent monthly data, HAC F is the more honest # # metric and preferred for monthly persistent data (NW lag = 4). # # ============================================================ # cat("--- 4. Instrument first-stage diagnostics ---\n\n") # Candidate instruments GOLD_INSTRUMENTS <- intersect( c("iv_gold_lbma_baseline", "iv_gold_lbma_core", "iv_gold_lbma_full"), names(lp_data) ) CES_INSTRUMENTS <- intersect( c("ces_abs", "ces_signed", "ces_neg_only", "d_ces_abs", "d_ces_signed", "d_ces_neg_only"), names(lp_data) ) GPR_INSTRUMENTS <- intersect( c("GPR", "gpr", "GPR_std", "gpr_std"), names(lp_data) ) SLOPE_INSTRUMENTS <- intersect( c("instrument_slope"), names(lp_data) ) # Optional financial-event gold IVs, if later created GOLD_FINANCIAL_INSTRUMENTS <- grep( "iv_gold.*financial|iv_gold_fin", names(lp_data), value = TRUE ) FS_SPECS <- list( ar_meu_innov = list( shocks = SHOCKS_MAIN, instruments = c(GOLD_INSTRUMENTS, CES_INSTRUMENTS, GPR_INSTRUMENTS) ), # CESIEUR is inside this PCA, so CESIEUR is NOT allowed as instrument. MEU_PCA_CES_innov = list( shocks = SHOCKS_PCA, instruments = c(GOLD_INSTRUMENTS, GPR_INSTRUMENTS) ), fu_uncertainty_innov = list( shocks = SHOCKS_MAIN, instruments = c(GOLD_INSTRUMENTS, GOLD_FINANCIAL_INSTRUMENTS, SLOPE_INSTRUMENTS) ), fu_riskpremium_innov = list( shocks = SHOCKS_MAIN, instruments = c(SLOPE_INSTRUMENTS, GOLD_INSTRUMENTS, GOLD_FINANCIAL_INSTRUMENTS) ) ) # Helper for lag blocks, local to first-stage section make_lag_block_fs <- function(data, vars, lags) { vars <- intersect(vars, names(data)) purrr::map_dfc(vars, function(v) { purrr::map_dfc(1:lags, function(L) { tibble(!!paste0(v, "_l", L) := dplyr::lag(data[[v]], L)) }) }) } run_fs_diag <- function(data, endog, instrument, shock_set, controls, outcome_for_lags = "ip_yoy", lags = N_LAGS, hac_lag = 4) { other_shocks <- setdiff(shock_set, endog) needed_raw <- c(endog, instrument, other_shocks, EXOG, outcome_for_lags, controls) missing_vars <- setdiff(needed_raw, names(data)) if (length(missing_vars) > 0) { stop("Missing variables in first-stage: ", paste(missing_vars, collapse = ", ")) } lag_y <- purrr::map_dfc(1:lags, function(L) { tibble(!!paste0(outcome_for_lags, "_l", L) := dplyr::lag(data[[outcome_for_lags]], L)) }) lag_controls <- make_lag_block_fs(data, controls, lags) reg_df <- bind_cols( data %>% select(all_of(c(endog, instrument, other_shocks, EXOG))), lag_y, lag_controls ) %>% drop_na() rhs <- paste( c(instrument, other_shocks, EXOG, names(lag_y), names(lag_controls)), collapse = " + " ) fit <- lm(as.formula(paste(endog, "~", rhs)), data = reg_df) s <- summary(fit) t_ols <- s$coefficients[instrument, "t value"] F_ols <- t_ols^2 vcov_hac <- NeweyWest(fit, lag = hac_lag, prewhite = FALSE, adjust = TRUE) ct_hac <- coeftest(fit, vcov. = vcov_hac) t_hac <- ct_hac[instrument, "t value"] F_hac <- t_hac^2 p_hac <- ct_hac[instrument, "Pr(>|t|)"] tibble( endog = endog, instrument = instrument, coef = s$coefficients[instrument, "Estimate"], t_ols = t_ols, F_ols = F_ols, t_hac = t_hac, F_hac = F_hac, p_hac = p_hac, n_obs = nobs(fit), verdict = case_when( F_hac >= 16 ~ "STRONG", F_hac >= 10 ~ "OK", F_hac >= 5 ~ "BORDERLINE", TRUE ~ "WEAK" ) ) } fs_results <- imap_dfr(FS_SPECS, function(spec, endog) { if (!endog %in% names(lp_data)) return(tibble()) instruments <- intersect(unique(spec$instruments), names(lp_data)) controls_fs <- get_controls("ip_yoy", spec = "baseline") map_dfr(instruments, function(iv) { tryCatch( run_fs_diag( data = lp_data, endog = endog, instrument = iv, shock_set = spec$shocks, controls = controls_fs, outcome_for_lags = "ip_yoy", lags = N_LAGS, hac_lag = 4 ), error = function(e) { cat("Skipped", endog, "~", iv, ":", conditionMessage(e), "\n") NULL } ) }) }) %>% mutate(across(where(is.numeric), ~ round(.x, 4))) %>% arrange(endog, desc(F_hac)) cat("First-stage diagnostics:\n") print(fs_results, n = Inf) write.xlsx( fs_results, file.path(TAB_MAIN, "first_stage_diagnostics.xlsx"), overwrite = TRUE ) strong_iv <- fs_results %>% filter(F_hac >= 10) if (nrow(strong_iv) > 0) { cat("\n*** Candidate instruments with HAC F >= 10 ***\n") print(strong_iv, n = Inf) cat(">>> Carry these only to the separate LP-IV script.\n\n") } else { cat("\nNo instrument passes HAC F >= 10.\n") cat("Proceeding with innovation-based LPs, not external-IV identification.\n\n") } save_table( fs_results %>% mutate(across(where(is.numeric), ~ round(.x, 3))), "first_stage_diagnostics", caption = "First-stage diagnostics: candidate instruments with HAC F >= 10", label = "tab:first_stage_diagnostics" ) # ============================================================ # # 5. HELPER FUNCTIONS # # ============================================================ # # 5a. Pre-create lag columns for a set of variables make_lag_block <- function(data, vars, lags) { purrr::map_dfc(vars, function(v) { purrr::map_dfc(1:lags, function(L) { tibble(!!paste0(v, "_l", L) := lag(data[[v]], L)) }) }) } # 5b. LP estimation # Methodology: # - All shocks enter contemporaneously (they are innovations) # - All state controls enter LAGGED only (predetermined) # - Outcome variable enters lagged (own dynamics) # - COVID dummy enters contemporaneously (exogenous) # - HAC standard errors: NW bandwidth = max(h, N_LAGS) # Following Jordà (2005): horizon-specific bandwidth # accounts for the MA(h) error structure induced by # multi-step projection # - 68% CI: ±1 SE (one-sigma band, common in macro LP) # - 90% CI: ±1.645 SE (standard significance band) run_lp <- function(data, outcome, shocks, controls, exo = NULL, H = 24, lags = N_LAGS) { n <- nrow(data) results <- vector("list", H + 1) # Pre-compute lag blocks (avoid lag() inside loop) lag_ctrls <- make_lag_block(data, controls, lags) for (h in 0:H) { # h-step-ahead outcome (Jordà 2005 projection) y_lead <- c(data[[outcome]][(h+1):n], rep(NA, h)) # Lags of outcome (own dynamics) lag_y <- purrr::map_dfc(1:lags, function(L) tibble(!!paste0(outcome,"_l",L) := lag(data[[outcome]], L))) # Contemporaneous exogenous dummies exo_df <- if (!is.null(exo) && length(exo) > 0) data[, exo, drop = FALSE] else tibble() reg_df <- bind_cols( tibble(y = y_lead), data[, shocks, drop = FALSE], # shocks at t (innovations) exo_df, # covid dummy at t lag_y, # outcome lags lag_ctrls # control lags ) %>% drop_na() rhs <- paste( c(shocks, names(exo_df), names(lag_y), names(lag_ctrls)), collapse = " + " ) fit <- lm(as.formula(paste("y ~", rhs)), data = reg_df) # Newey-West: bandwidth = max(h, N_LAGS) # Rationale: LP at horizon h has MA(h) errors by construction, # requiring at least h lags in the HAC estimator nw_bw <- max(h, lags) vcov_nw <- NeweyWest(fit, lag = nw_bw, prewhite = FALSE, adjust = TRUE) se_nw <- sqrt(diag(vcov_nw)) cf <- coef(fit) results[[h+1]] <- purrr::map_dfr(shocks, function(s) { tibble( horizon = h, shock = s, coef = cf[s], se = se_nw[s], tstat = cf[s] / se_nw[s], pval = 2 * pt(-abs(cf[s]/se_nw[s]), df = nobs(fit) - length(cf)), ci90_lo = cf[s] - 1.645 * se_nw[s], ci90_hi = cf[s] + 1.645 * se_nw[s], ci68_lo = cf[s] - 1.000 * se_nw[s], ci68_hi = cf[s] + 1.000 * se_nw[s], sig90 = abs(cf[s]/se_nw[s]) > 1.645, sig68 = abs(cf[s]/se_nw[s]) > 1.000, n_obs = nobs(fit), adj_r2 = summary(fit)$adj.r.squared, nw_bw = nw_bw ) }) } bind_rows(results) } # 5c. Coefficient table at key horizons print_irf_table <- function(irf_data, label = "") { cat("\n=== IRF coefficients:", label, "===\n") cat(sprintf("%-6s %-26s %8s %8s %8s %6s\n", "h", "Shock", "Coef", "SE", "t-stat", "Sig")) cat(strrep("-", 68), "\n") irf_data %>% filter(horizon %in% c(0, 3, 6, 9, 12, 18, 24)) %>% arrange(horizon, shock) %>% mutate(sig = case_when(pval<.01~"***", pval<.05~"**", pval<.10~"*", TRUE~"")) %>% pwalk(function(horizon, shock, coef, se, tstat, sig, ...) { cat(sprintf("h=%-4d %-26s %8.4f %8.4f %8.3f %6s\n", horizon, shock, coef, se, tstat, sig)) }) cat("Significance: *** p<0.01 ** p<0.05 * p<0.10\n") } # 5d. Single IRF panel plot_irf <- function(irf_data, shock_name, title, subtitle = NULL, color = "#2166ac", y_label = "Response (pp)") { df <- irf_data %>% filter(shock == shock_name) ggplot(df, aes(x = horizon)) + geom_hline(yintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.5) + geom_ribbon(aes(ymin = ci90_lo, ymax = ci90_hi), fill = color, alpha = 0.18) + geom_ribbon(aes(ymin = ci68_lo, ymax = ci68_hi), fill = color, alpha = 0.30) + geom_line(aes(y = coef), color = color, linewidth = 1.0) + geom_point(aes(y = coef, shape = sig90, size = sig90), color = color) + scale_shape_manual(values = c("FALSE"=1, "TRUE"=19)) + scale_size_manual( values = c("FALSE"=1.5, "TRUE"=2.5)) + scale_x_continuous(breaks = seq(0, 24, by = 4)) + labs(title = title, subtitle = subtitle, x = "Horizon (months)", y = y_label) + theme_minimal() + theme(legend.position = "none", plot.title = element_text(size=10, face="bold"), plot.subtitle = element_text(size=8, color="grey40"), panel.grid.minor = element_blank()) } # 5e. Three-shock panel for one outcome plot_three_shocks <- function(irf_data, outcome_label, y_label, shocks = SHOCKS_MAIN) { plots <- map(shocks, function(s) { plot_irf( irf_data, shock_name = s, title = paste(shock_labels[s], "\u2192", outcome_label), subtitle = "Newey-West HAC. Shaded: 68% and 90% CI. \u25CF = sig. at 90%.", color = shock_colors[s], y_label = y_label ) }) wrap_plots(plots, ncol = 1) } # 5f. MU shock comparison: main vs PCA robustness plot_mu_comparison <- function(irf_main, irf_pca, outcome_label, y_label) { d_main <- irf_main %>% filter(shock == "ar_meu_innov") %>% mutate(spec = "Main (MEU innov.)") d_pca <- irf_pca %>% filter(shock == "MEU_PCA_CES_innov") %>% mutate(spec = "Robustness (PCA-CES innov.)") bind_rows(d_main, d_pca) %>% ggplot(aes(x = horizon)) + geom_hline(yintercept = 0, linetype = "dashed", color = "grey50", linewidth = 0.5) + geom_ribbon(aes(ymin = ci90_lo, ymax = ci90_hi, fill = spec), alpha = 0.18) + geom_line(aes(y = coef, color = spec), linewidth = 1.0) + geom_point(aes(y = coef, color = spec, shape = sig90, size = sig90)) + scale_shape_manual(values = c("FALSE"=1, "TRUE"=19)) + scale_size_manual( values = c("FALSE"=1.5,"TRUE"=2.5)) + scale_color_manual(values = c( "Main (MEU innov.)" = "#2166ac", "Robustness (PCA-CES innov.)" = "#4393c3" )) + scale_fill_manual(values = c( "Main (MEU innov.)" = "#2166ac", "Robustness (PCA-CES innov.)" = "#4393c3" )) + scale_x_continuous(breaks = seq(0, 24, by = 4)) + labs( title = paste("MU shock \u2192", outcome_label), subtitle = "Blue: main spec. Light blue: PCA-CES robustness. 90% CI.", x = "Horizon (months)", y = y_label, color = NULL, fill = NULL ) + theme_minimal() + theme(legend.position = "bottom", plot.title = element_text(size=10, face="bold"), plot.subtitle = element_text(size=8, color="grey40"), panel.grid.minor = element_blank()) } # ============================================================ # # 6. RUN LP — MAIN SPECIFICATION # # Outcomes: # # Primary: ip_yoy, inflation, unemployment # # Secondary: capital_goods_growth, nfc_growth # # # # Note on outcome choice: # # ip_yoy: standard real activity measure in LP uncertainty # # d_unemployment: first difference — same rationale + # # Okun's law links it directly to output gap # # capital_goods_growth: investment channel — directly tests # # the wait-and-see mechanism (Bloom 2009) # # nfc_growth: credit channel — BIS/ECB research on # # uncertainty → corporate credit contraction # # ============================================================ # # ============================================================ # # 6. RUN LPs — MAIN AND PCA_CES ROBUSTNESS # # ============================================================ # cat("\n\n=== 6. RUNNING BASELINE LPs ===\n") outcome_labels <- c( ip_yoy = "IP Growth (YoY %)", inflation = "Inflation rate", unemployment = "Unemployment rate", capital_goods_growth = "Capital Goods Growth", nfc_growth = "NFC Loan Growth", ciss_std = "CISS financial stress", ciss = "CISS financial stress", credit_spreads_std = "Credit spreads", credit_spreads = "Credit spreads", d_credit_spreads = "Δ Credit spreads" ) outcome_ylabels <- c( ip_yoy = "Response (pp)", inflation = "Response (pp)", unemployment = "Response (pp)", capital_goods_growth = "Response (pp)", nfc_growth = "Response (pp)", ciss_std = "Response (std units)", ciss = "Response", credit_spreads_std = "Response (std units)", credit_spreads = "Response (pp/spread units)", d_credit_spreads = "Response (change)" ) run_outcome_suite <- function(data, shocks, shock_spec_label, control_spec = "baseline", lags = N_LAGS, verbose = TRUE) { res <- list() controls_used <- list() for (y in OUTCOMES) { controls_y <- get_controls(y, spec = control_spec) controls_used[[y]] <- controls_y if (verbose) { cat("\n======================================================\n") cat("Shock specification:", shock_spec_label, "\n") cat("Outcome:", y, "-", outcome_labels[y], "\n") cat("Control spec:", control_spec, "\n") cat("Controls used:", paste(controls_y, collapse = ", "), "\n") cat("Sample before run_lp:", nrow(data), "rows\n") cat("LP lags:", lags, "\n") cat("======================================================\n") } res[[y]] <- run_lp( data = data, outcome = y, shocks = shocks, controls = controls_y, exo = EXOG, H = H, lags = lags ) %>% mutate( outcome = y, outcome_label = outcome_labels[y], shock_spec = shock_spec_label, control_spec = control_spec, controls_used = paste(controls_y, collapse = " + ") ) if (verbose) { print_irf_table( res[[y]], paste0(outcome_labels[y], " [", shock_spec_label, "]") ) } } attr(res, "controls_used") <- controls_used return(res) } # Main specification: MEU + BH FU innovations irfs_main <- run_outcome_suite( data = lp_data, shocks = SHOCKS_MAIN, shock_spec_label = "main_meu", control_spec = "baseline", lags = N_LAGS, verbose = TRUE ) # PCA-CES robustness: replace MU shock, FU shocks unchanged lp_pca <- lp_data %>% filter(!is.na(MEU_PCA_CES_innov)) irfs_pca <- run_outcome_suite( data = lp_pca, shocks = SHOCKS_PCA, shock_spec_label = paste0("pca_ces_L", N_LAGS), control_spec = "baseline", lags = N_LAGS, verbose = TRUE ) # Optional extended controls with ESI RUN_EXTENDED_CONTROLS <- length(ESI_CONTROL) > 0 if (RUN_EXTENDED_CONTROLS) { cat("\n\n=== 6B. RUNNING EXTENDED-CONTROL ROBUSTNESS LPs ===\n") cat("Extended control added:", ESI_CONTROL[1], "\n") irfs_main_ext <- run_outcome_suite( data = lp_data, shocks = SHOCKS_MAIN, shock_spec_label = "main_meu_extended", control_spec = "extended" ) irfs_pca_ext <- run_outcome_suite( data = lp_pca, shocks = SHOCKS_PCA, shock_spec_label = "pca_ces_extended", control_spec = "extended" ) } # ============================================================ # # ANNEX: Contemporaneous shock contribution to IP growth (h=0) # # ============================================================ # cat("\n--- Annex: Contemporaneous shock decomposition (IP, h=0) ---\n") # Extract h=0 coefficients for IP growth from irfs_main ip_h0 <- irfs_main %>% bind_rows() %>% filter(outcome == "ip_yoy", horizon == 0) beta_meu <- ip_h0 %>% filter(shock == "ar_meu_innov") %>% pull(coef) beta_fu <- ip_h0 %>% filter(shock == "fu_uncertainty_innov") %>% pull(coef) beta_vrp <- ip_h0 %>% filter(shock == "fu_riskpremium_innov") %>% pull(coef) decomp_df <- lp_data %>% select(date, ar_meu_innov, fu_uncertainty_innov, fu_riskpremium_innov) %>% drop_na() %>% transmute( date, MEU = beta_meu * ar_meu_innov, FU = beta_fu * fu_uncertainty_innov, VRP = beta_vrp * fu_riskpremium_innov ) %>% pivot_longer(-date, names_to = "shock", values_to = "contribution") %>% mutate(shock = factor(shock, levels = c("MEU", "FU", "VRP"))) p_decomp <- ggplot(decomp_df, aes(x = date, y = contribution, fill = shock)) + geom_col(position = "stack", width = 31) + geom_hline(yintercept = 0, color = "grey30", linewidth = 0.4) + scale_fill_manual(values = c(MEU = "#2166ac", FU = "#d6604d", VRP = "#4dac26")) + labs( title = "Contemporaneous Shock Contributions to IP Growth (h = 0)", subtitle = "beta_h0 x innovation, by uncertainty component. Illustrative decomposition, not a structural historical decomposition.", x = NULL, y = "Contribution (pp)", fill = NULL ) + theme_minimal() + theme( plot.title = element_text(face = "bold", size = 11), plot.subtitle = element_text(size = 8, color = "grey35"), legend.position = "bottom" ) save_fig(p_decomp, "annex_ip_contemporaneous_decomp.pdf", dir = FIG_MAIN, width = 11, height = 5) cat("Saved: annex_ip_contemporaneous_decomp.pdf\n") # ============================================================ # # 6C. LAG-LENGTH ROBUSTNESS — PRIMARY OUTCOME ONLY # # ============================================================ # cat("\n\n=== 6C. LAG-LENGTH ROBUSTNESS: IP_YOY ===\n") lag_robust_ip <- purrr::map_dfr(LAG_ROBUSTNESS_SET, function(L) { cat("\nRunning IP LP with", L, "lags...\n") run_lp( data = lp_data, outcome = "ip_yoy", shocks = SHOCKS_MAIN, controls = get_controls("ip_yoy", spec = "baseline"), exo = EXOG, H = H, lags = L ) %>% mutate(lags = L) }) save_table( lag_robust_ip %>% mutate(across(where(is.numeric), ~ round(.x, 3))), "lag_robustness_ip", caption = " Lag Robustness IP", label = "tab:lag_robustness_ip" ) # Quick console comparison for key horizons cat("\nLag robustness, selected horizons:\n") lag_robust_ip %>% filter(horizon %in% c(0, 3, 6, 9, 12, 18, 24)) %>% select(lags, horizon, shock, coef, se, pval) %>% arrange(shock, horizon, lags) %>% mutate(across(c(coef, se, pval), ~ round(.x, 4))) %>% print(n = Inf) p_lag_robust_mu <- lag_robust_ip %>% filter(shock == "ar_meu_innov") %>% ggplot(aes(x = horizon, y = coef, color = factor(lags))) + geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") + geom_line(linewidth = 1) + geom_point(size = 1.5) + labs( title = "Lag-length robustness: MU shock → IP growth", subtitle = "Comparison across selected LP lag lengths", x = "Horizon (months)", y = "Response (pp)", color = "LP lags" ) + theme_minimal() + theme(legend.position = "bottom") save_fig( p_lag_robust_mu, "lag_robustness_mu_ip.pdf", FIG_MAIN, width = 8, height = 5 ) # ============================================================ # # 7. SAVE RESULTS TO TABLES # # ============================================================ # cat("\n--- 8. Saving IRF tables ---\n") save_irf_list <- function(irf_list, dir_out, prefix) { for (nm in names(irf_list)) { write.xlsx( irf_list[[nm]], file.path(dir_out, paste0(prefix, "_", nm, ".xlsx")), overwrite = TRUE ) } } save_irf_list(irfs_main, TAB_MAIN, "irf_main") save_irf_list(irfs_pca, TAB_PCA, "irf_pca") if (exists("irfs_main_ext")) { save_irf_list(irfs_main_ext, TAB_MAIN, "irf_main_extended") } if (exists("irfs_pca_ext")) { save_irf_list(irfs_pca_ext, TAB_PCA, "irf_pca_extended") } # Keep old object names for plotting compatibility irf_ip <- irfs_main[["ip_yoy"]] irf_inf <- irfs_main[["inflation"]] irf_unemp <- irfs_main[["unemployment"]] irf_capg <- irfs_main[["capital_goods_growth"]] irf_nfc <- irfs_main[["nfc_growth"]] irf_ciss <- irfs_main[[CISS_OUTCOME]] irf_credit <- irfs_main[[CREDIT_OUTCOME]] irf_ip_pca <- irfs_pca[["ip_yoy"]] irf_inf_pca <- irfs_pca[["inflation"]] irf_unemp_pca <- irfs_pca[["unemployment"]] irf_ciss_pca <- irfs_pca[[CISS_OUTCOME]] irf_credit_pca <- irfs_pca[[CREDIT_OUTCOME]] cat("Tables saved.\n") # ============================================================ # # 9. IRF PLOTS — MAIN SPECIFICATION # # ============================================================ # cat("\n--- 9. Plotting main IRFs ---\n") # 9.1 All shocks × IP growth p_ip <- plot_three_shocks(irf_ip, "IP Growth (YoY %)", "Response (pp)") save_fig(p_ip, "irf_main_ip.pdf", FIG_MAIN, width=10, height=15) # 9.2 All shocks × Inflation p_inf <- plot_three_shocks(irf_inf, "Inflation (MoM change)", "Response (pp)") save_fig(p_inf, "irf_main_inflation.pdf", FIG_MAIN, width=10, height=15) # 9.3 All shocks × Unemployment p_unemp <- plot_three_shocks(irf_unemp, "\u0394 Unemployment", "Response (pp change)") save_fig(p_unemp, "irf_main_unemployment.pdf", FIG_MAIN, width=10, height=15) # 9.4 All shocks × Capital Goods p_capg <- plot_three_shocks(irf_capg, "Capital Goods Growth", "Response (pp)") save_fig(p_capg, "irf_main_capital_goods.pdf", FIG_MAIN, width=10, height=15) # 9.5 All shocks × NFC loans p_nfc <- plot_three_shocks(irf_nfc, "NFC Loan Growth", "Response (pp)") save_fig(p_nfc, "irf_main_nfc.pdf", FIG_MAIN, width=10, height=15) # 9.6 CISS financial stress p_ciss <- plot_three_shocks( irf_ciss, outcome_labels[CISS_OUTCOME], outcome_ylabels[CISS_OUTCOME] ) save_fig(p_ciss, "irf_main_ciss.pdf", FIG_MAIN, width = 10, height = 15) # 9.7 Credit spreads p_credit <- plot_three_shocks( irf_credit, outcome_labels[CREDIT_OUTCOME], outcome_ylabels[CREDIT_OUTCOME] ) save_fig(p_credit, "irf_main_credit_spreads.pdf", FIG_MAIN, width = 10, height = 15) # 9.8 Three-outcome summary for each shock (thesis-ready panels) # For each shock: IP / Inflation / Unemployment side by side for (s in SHOCKS_MAIN) { p_summary <- ( plot_irf(irf_ip, s, paste(shock_labels[s], "\u2192 IP Growth"), color = shock_colors[s], y_label = "Response (pp)") + plot_irf(irf_inf, s, paste(shock_labels[s], "\u2192 Inflation"), color = shock_colors[s], y_label = "Response (pp)") + plot_irf(irf_unemp, s, paste(shock_labels[s], "\u2192 Unemployment"), color = shock_colors[s], y_label = "Response (pp change)") ) + plot_annotation( title = "LP Impulse Responses — Main Specification", subtitle = "BH FU decomposition. Innovations as shocks. NW-HAC. 68% + 90% CI." ) fname <- paste0("irf_summary_", gsub("_innov","",s), ".pdf") save_fig(p_summary, fname, FIG_MAIN, width=15, height=5) } # ============================================================ # # 11. SPECIFICATION SUMMARY # # ============================================================ # cat("\n=== SPECIFICATION SUMMARY ===\n") cat("Main specification:\n") cat(" MU shock: ar_meu_innov\n") cat(" FU level: fu_uncertainty_innov (BH expected physical variance)\n") cat(" FU prem.: fu_riskpremium_innov (BH variance risk premium)\n") cat(" Identification: innovation-based LP, not external-IV identification.\n") cat(" Shocks enter contemporaneously; controls enter lagged.\n\n") cat("PCA-CES robustness:\n") cat(" MU shock: MEU_PCA_CES_innov\n") cat(" Meaning: PCA(MEU + EPU + CESIEUR, excluding GPR)\n") cat(" Note: CESIEUR-augmented macro uncertainty robustness factor.\n") cat(" FU shocks: unchanged.\n\n") cat("Controls:\n") cat(" Baseline controls are outcome-specific.\n") cat(" Inflation/unemployment/capital-goods equations include lagged IP.\n") cat(" CISS and credit-spread equations remove the corresponding financial control.\n") cat(" ESI extended robustness: ", ifelse(RUN_EXTENDED_CONTROLS, ESI_CONTROL[1], "not run / not available"), "\n\n") cat("Outcomes:\n") cat(" ", paste(OUTCOMES, collapse = ", "), "\n\n") cat("LP settings:\n") cat(" Lags: ", N_LAGS, "\n") cat(" Horizons: 0 to", H, "months\n") cat(" SE: Newey-West HAC, bandwidth = max(h, N_LAGS)\n") cat(" CI: 68% and 90% bands\n\n") cat("Output:\n") cat(" Figures: ", FIG_MAIN, "/\n") cat(" Figures: ", FIG_PCA, "/\n") cat(" Tables: ", TAB_MAIN, "/\n") cat(" Tables: ", TAB_PCA, "/\n") cat("\n=== I_LP_BASELINE COMPLETE ===\n") ## EXTRA CHECK TO SEE WHY IT CHANGES ## # Same-sample comparison common_sample <- lp_data %>% filter( !is.na(ar_meu_innov), !is.na(MEU_PCA_CES_innov), !is.na(fu_uncertainty_innov), !is.na(fu_riskpremium_innov) ) irf_main_common <- run_lp( data = common_sample, outcome = "ip_yoy", shocks = SHOCKS_MAIN, controls = get_controls("ip_yoy", spec = "baseline"), exo = EXOG, H = H, lags = N_LAGS ) irf_pca_common <- run_lp( data = common_sample, outcome = "ip_yoy", shocks = SHOCKS_PCA, controls = get_controls("ip_yoy", spec = "baseline"), exo = EXOG, H = H, lags = N_LAGS ) common_sample %>% select( ar_meu_innov, MEU_PCA_CES_innov, fu_uncertainty_innov, fu_riskpremium_innov ) %>% cor(use = "complete.obs") %>% round(3) save_table( as.data.frame(round(cor_mat, 3)), "shock_correlations", caption = "Pairwise correlations between uncertainty innovations", label = "tab:shock_correlations" ) cat("\nInterpretation notes:\n") cat(" MEU-FU correlation 0.414: moderate overlap, not collinear\n") cat(" MEU-VRP correlation 0.043: near-orthogonal, strong separation\n") cat(" FU-VRP correlation -0.259: BH decomposition working as expected\n") cat(" All below 0.45: joint LP entry is well-identified\n\n") # ============================================================ # # SECTION 12: THESIS-READY DASHBOARD FIGURES # # ============================================================ # # ── Helper: single IRF panel, paper style ───────────────────── plot_irf_clean <- function(irf_data, shock_name, outcome_label, color, y_label = "Response (pp)") { df <- irf_data %>% filter(shock == shock_name) ggplot(df, aes(x = horizon)) + geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3, color = "grey40") + geom_ribbon(aes(ymin = ci90_lo, ymax = ci90_hi), fill = color, alpha = 0.15) + geom_ribbon(aes(ymin = ci68_lo, ymax = ci68_hi), fill = color, alpha = 0.28) + geom_line(aes(y = coef), color = color, linewidth = 0.8) + geom_point(aes(y = coef, shape = sig90, size = sig90), color = color) + scale_shape_manual(values = c("FALSE" = 1, "TRUE" = 19)) + scale_size_manual(values = c("FALSE" = 1.2, "TRUE" = 2.2)) + scale_x_continuous(breaks = seq(0, 24, by = 6)) + labs(x = "Horizon (months)", y = y_label, title = outcome_label) + theme_classic(base_size = 9) + theme( legend.position = "none", plot.title = element_text(size = 9, face = "bold"), axis.title = element_text(size = 8), axis.text = element_text(size = 7.5), panel.grid.major.y = element_line(linewidth = 0.2, color = "grey90"), plot.margin = margin(4, 6, 4, 4) ) } # ── Helper: build one shock-row of 3 panels ─────────────────── build_shock_row <- function(irf_list, shock, outcomes, labels, ylabels, color) { purrr::map(outcomes, function(y) { plot_irf_clean( irf_data = irf_list[[y]], shock_name = shock, outcome_label = labels[y], color = color, y_label = ylabels[y] ) }) } # ── Outcome sets ────────────────────────────────────────────── MACRO_OUTCOMES <- c("ip_yoy", "inflation", "unemployment") FINANCIAL_OUTCOMES <- c( "capital_goods_growth", "nfc_growth", CISS_OUTCOME, CREDIT_OUTCOME ) %>% purrr::discard(is.na) %>% intersect(names(lp_data)) # Use first 3 financial outcomes for the 3x3 grid FINANCIAL_OUTCOMES_3 <- FINANCIAL_OUTCOMES[1:min(3, length(FINANCIAL_OUTCOMES))] short_labels <- c( ip_yoy = "IP Growth", inflation = "Inflation", unemployment = "Unemployment", capital_goods_growth = "Capital Goods", nfc_growth = "NFC Loans", ciss_std = "CISS Stress", ciss = "CISS Stress", credit_spreads_std = "Credit Spreads", credit_spreads = "Credit Spreads", d_credit_spreads = "\u0394 Credit Spreads" ) short_ylabels <- c( ip_yoy = "Response (pp)", inflation = "Response (pp)", unemployment = "Response (pp)", capital_goods_growth = "Response (pp)", nfc_growth = "Response (pp)", ciss_std = "Response (std)", ciss = "Response", credit_spreads_std = "Response (std)", credit_spreads = "Response", d_credit_spreads = "Response" ) shocks_ordered <- c("ar_meu_innov", "fu_uncertainty_innov", "fu_riskpremium_innov") shock_row_labels <- c( ar_meu_innov = "MU shock", fu_uncertainty_innov = "FU shock", fu_riskpremium_innov = "VRP shock" ) # ── Function: build 3x3 dashboard ──────────────────────────── build_3x3_dashboard <- function(irf_list, outcomes, title_text, filename, fig_dir = FIG_MAIN) { all_panels <- list() for (s in shocks_ordered) { row_panels <- build_shock_row( irf_list = irf_list, shock = s, outcomes = outcomes, labels = short_labels, ylabels = short_ylabels, color = shock_colors[s] ) all_panels <- c(all_panels, row_panels) } # Row labels as left-side annotations row_label_grobs <- purrr::map(shocks_ordered, function(s) { ggplot() + annotate("text", x = 0.5, y = 0.5, label = shock_row_labels[s], angle = 90, fontface = "bold", size = 3.2, color = shock_colors[s]) + theme_void() + theme(plot.margin = margin(0, 2, 0, 2)) }) # Interleave: label | p1 | p2 | p3 for each row full_layout <- list() for (i in seq_along(shocks_ordered)) { full_layout <- c( full_layout, list(row_label_grobs[[i]]), all_panels[((i-1)*length(outcomes)+1):(i*length(outcomes))] ) } n_cols <- length(outcomes) + 1 # +1 for row label n_rows <- length(shocks_ordered) col_widths <- c(0.12, rep(1, length(outcomes))) p_final <- patchwork::wrap_plots(full_layout, ncol = n_cols, widths = col_widths) + patchwork::plot_annotation( title = title_text, subtitle = "Baseline LP. NW-HAC. Shaded: 68% (dark) and 90% (light) CI. \u25CF = sig. at 90%.", theme = theme( plot.title = element_text(face = "bold", size = 11), plot.subtitle = element_text(size = 8.5, color = "grey35") ) ) save_fig(p_final, filename, dir = fig_dir, width = 11, height = 9) invisible(p_final) } # ============================================================ # # 12A. MAIN TEXT: 3x3 MACRO DASHBOARD # # ============================================================ # cat("\n--- 12A. Building macro 3x3 dashboard ---\n") build_3x3_dashboard( irf_list = irfs_main, outcomes = MACRO_OUTCOMES, title_text = "Baseline LP: Three Uncertainty Shocks \u00D7 Macroeconomic Outcomes", filename = "dashboard_3x3_macro.pdf" ) # ============================================================ # # 12B. MAIN TEXT: 3x3 FINANCIAL DASHBOARD # # ============================================================ # cat("\n--- 12B. Building financial 3x3 dashboard ---\n") if (length(FINANCIAL_OUTCOMES_3) == 3) { build_3x3_dashboard( irf_list = irfs_main, outcomes = FINANCIAL_OUTCOMES_3, title_text = "Baseline LP: Three Uncertainty Shocks \u00D7 Financial Outcomes", filename = "dashboard_3x3_financial.pdf" ) } else { cat("Only", length(FINANCIAL_OUTCOMES_3), "financial outcomes available — need 3 for dashboard.\n") cat("Available:", paste(FINANCIAL_OUTCOMES_3, collapse = ", "), "\n") } # ============================================================ # # 12C. ANNEX: MAIN vs PCA-CES COMPARISON # # One figure per shock. Each figure: 3 macro outcomes, # # main spec (solid) vs PCA-CES robustness (dashed). # # Uses same colour scheme as existing plot_mu_comparison. # # ============================================================ # cat("\n--- 12C. Annex: main vs PCA-CES comparison ---\n") # Colour pairs per shock: main = full colour, robust = lighter shade pca_color_pairs <- list( ar_meu_innov = c(main = "#2166ac", robust = "#6baed6"), fu_uncertainty_innov = c(main = "#d6604d", robust = "#f4a582"), fu_riskpremium_innov = c(main = "#4dac26", robust = "#a1d76a") ) # PCA shock name lookup pca_shock_name <- function(s) { if (s == "ar_meu_innov") "MEU_PCA_CES_innov" else s } build_pca_comparison_panel <- function(s, outcomes, irf_main_list, irf_pca_list) { col_pairs <- list( ar_meu_innov = c(main = "#2166ac", robust = "#6baed6"), fu_uncertainty_innov = c(main = "#d6604d", robust = "#f4a582"), fu_riskpremium_innov = c(main = "#4dac26", robust = "#a1d76a") ) col_main <- col_pairs[[s]]["main"] col_pca <- col_pairs[[s]]["robust"] pca_s <- if (s == "ar_meu_innov") "MEU_PCA_CES_innov" else s panels <- purrr::map(outcomes, function(y) { d_main <- irf_main_list[[y]] %>% filter(shock == s) %>% mutate(spec = "Main") d_pca <- irf_pca_list[[y]] %>% filter(shock == pca_s) %>% mutate(spec = "PCA-CES") if (nrow(d_main) == 0 || nrow(d_pca) == 0) return(NULL) ggplot() + # PCA ribbon + line — drawn first (underneath) geom_ribbon(data = d_pca, aes(x = horizon, ymin = ci90_lo, ymax = ci90_hi), fill = as.character(col_pca), alpha = 0.20) + geom_line(data = d_pca, aes(x = horizon, y = coef), color = as.character(col_pca), linetype = "dashed", linewidth = 0.9) + # Main ribbon + line — drawn on top geom_ribbon(data = d_main, aes(x = horizon, ymin = ci90_lo, ymax = ci90_hi), fill = as.character(col_main), alpha = 0.20) + geom_ribbon(data = d_main, aes(x = horizon, ymin = ci68_lo, ymax = ci68_hi), fill = as.character(col_main), alpha = 0.28) + geom_line(data = d_main, aes(x = horizon, y = coef), color = as.character(col_main), linetype = "solid", linewidth = 1.0) + geom_point(data = d_main %>% filter(sig90 == TRUE), aes(x = horizon, y = coef), color = as.character(col_main), shape = 19, size = 2.2) + geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3, color = "grey40") + scale_x_continuous(breaks = seq(0, 24, by = 4)) + labs( title = paste0(shock_row_labels[s], " \u2192 ", short_labels[y]), x = "Horizon (months)", y = short_ylabels[y] ) + theme_minimal(base_size = 9.5) + theme( legend.position = "none", plot.title = element_text(size = 10, face = "bold"), axis.title = element_text(size = 8.5), axis.text = element_text(size = 8), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(linewidth = 0.25, color = "grey88"), plot.margin = margin(4, 6, 4, 4) ) }) %>% purrr::compact() # Shared legend as annotation text instead of ggplot legend patchwork::wrap_plots(panels, ncol = 3) + patchwork::plot_annotation( title = paste0(shock_row_labels[s], " \u2014 Main vs PCA-CES robustness"), subtitle = paste0( "Solid (", col_main["main"], "): main specification. ", "Dashed (lighter): PCA-CES robustness. ", "Dark shading: 68% CI. Light shading: 90% CI." ), theme = theme( plot.title = element_text(face = "bold", size = 11), plot.subtitle = element_text(size = 8.5, color = "grey35") ) ) } # Run 3 annex figures (one per shock, 3 macro outcomes each) for (s in shocks_ordered) { p_annex <- build_pca_comparison_panel( s = s, outcomes = MACRO_OUTCOMES, irf_main_list = irfs_main, irf_pca_list = irfs_pca ) fname <- paste0("annex_pca_comparison_", gsub("_innov", "", s), "_macro.pdf") save_fig(p_annex, fname, dir = FIG_PCA, width = 13, height = 5) cat("Saved:", fname, "\n") } # ============================================================ # # 12D. ANNEX: ESI ROBUSTNESS — 3 shocks x 3 macro outcomes # # Adds ESI as additional control. Tests whether including # # consumer/business sentiment changes the baseline IRFs. # # ESI is NOT treated as an outcome here — it is a control # # robustness check only. # # ============================================================ # cat("\n--- 12D. ESI robustness ---\n") if (length(ESI_CONTROL) > 0) { cat("ESI found:", ESI_CONTROL[1], "\n") irfs_esi <- purrr::map(MACRO_OUTCOMES, function(y) { controls_esi <- unique(c( get_controls(y, spec = "baseline"), ESI_CONTROL[1] )) controls_esi <- intersect(controls_esi, names(lp_data)) run_lp( data = lp_data, outcome = y, shocks = SHOCKS_MAIN, controls = controls_esi, exo = EXOG, H = H, lags = N_LAGS ) }) names(irfs_esi) <- MACRO_OUTCOMES esi_color_pairs <- list( ar_meu_innov = c(base = "#2166ac", esi = "#6baed6"), fu_uncertainty_innov = c(base = "#d6604d", esi = "#f4a582"), fu_riskpremium_innov = c(base = "#4dac26", esi = "#a1d76a") ) esi_panels <- list() for (s in shocks_ordered) { for (y in MACRO_OUTCOMES) { col_base <- as.character(esi_color_pairs[[s]]["base"]) col_esi <- as.character(esi_color_pairs[[s]]["esi"]) d_base <- irfs_main[[y]] %>% filter(shock == s) d_esi <- irfs_esi[[y]] %>% filter(shock == s) esi_panels[[paste(s, y, sep = "_")]] <- ggplot() + geom_ribbon(data = d_esi, aes(x = horizon, ymin = ci90_lo, ymax = ci90_hi), fill = col_esi, alpha = 0.20) + geom_line(data = d_esi, aes(x = horizon, y = coef), color = col_esi, linetype = "dashed", linewidth = 0.9) + geom_ribbon(data = d_base, aes(x = horizon, ymin = ci90_lo, ymax = ci90_hi), fill = col_base, alpha = 0.20) + geom_ribbon(data = d_base, aes(x = horizon, ymin = ci68_lo, ymax = ci68_hi), fill = col_base, alpha = 0.28) + geom_line(data = d_base, aes(x = horizon, y = coef), color = col_base, linetype = "solid", linewidth = 1.0) + geom_point(data = d_base %>% filter(sig90 == TRUE), aes(x = horizon, y = coef), color = col_base, shape = 19, size = 2.2) + geom_hline(yintercept = 0, linetype = "dashed", linewidth = 0.3, color = "grey40") + scale_x_continuous(breaks = seq(0, 24, by = 4)) + labs( title = paste0(shock_row_labels[s], " \u2192 ", short_labels[y]), x = "Horizon (months)", y = short_ylabels[y] ) + theme_minimal(base_size = 9.5) + theme( legend.position = "none", plot.title = element_text(size = 10, face = "bold"), axis.title = element_text(size = 8.5), axis.text = element_text(size = 8), panel.grid.minor = element_blank(), panel.grid.major.x = element_blank(), panel.grid.major.y = element_line(linewidth = 0.25, color = "grey88"), plot.margin = margin(4, 6, 4, 4) ) } } p_esi <- patchwork::wrap_plots(esi_panels, ncol = 3) + patchwork::plot_annotation( title = "ESI robustness: baseline vs ESI-augmented controls", subtitle = paste0( "Solid darker: baseline specification. ", "Dashed lighter: baseline + ", ESI_CONTROL[1], " as additional control. 90% CI shaded." ), theme = theme( plot.title = element_text(face = "bold", size = 11), plot.subtitle = element_text(size = 8.5, color = "grey35") ) ) save_fig(p_esi, "annex_esi_robustness.pdf", dir = FIG_PCA, width = 13, height = 12) cat("ESI robustness saved.\n") } else { cat("ESI not found. Skipping.\n") } cat("\n=== SECTION 12 COMPLETE ===\n") cat("Main: dashboard_3x3_macro.pdf\n") cat("Main: dashboard_3x3_financial.pdf\n") cat("Annex: annex_pca_comparison_ar_meu_macro.pdf\n") cat("Annex: annex_pca_comparison_fu_uncertainty_macro.pdf\n") cat("Annex: annex_pca_comparison_fu_riskpremium_macro.pdf\n") cat("Annex: annex_esi_robustness.pdf\n") # ============================================================ # # APPENDIX TABLES: IRF COEFFICIENTS AND FULL IRF GRIDS # # ============================================================ # cat("\n--- Creating LaTeX-ready IRF appendix tables ---\n") TEX_MAIN <- file.path(TAB_MAIN, "latex") dir.create(TEX_MAIN, recursive = TRUE, showWarnings = FALSE) if (!requireNamespace("knitr", quietly = TRUE)) install.packages("knitr") library(knitr) # ------------------------------------------------------------ # # Helper functions # ------------------------------------------------------------ # star_fun <- function(p) { case_when( p < 0.01 ~ "***", p < 0.05 ~ "**", p < 0.10 ~ "*", TRUE ~ "" ) } fmt_coef <- function(coef, pval, digits = 3) { paste0(sprintf(paste0("%.", digits, "f"), coef), star_fun(pval)) } fmt_coef_se <- function(coef, se, pval, digits = 3) { paste0( sprintf(paste0("%.", digits, "f"), coef), star_fun(pval), " (", sprintf(paste0("%.", digits, "f"), se), ")" ) } write_tex <- function(tab, filename, caption, label, digits = 3) { tex <- knitr::kable( tab, format = "latex", booktabs = TRUE, caption = caption, label = label, digits = digits, escape = FALSE ) writeLines(tex, file.path(TEX_MAIN, paste0(filename, ".tex"))) cat("Saved:", file.path(TEX_MAIN, paste0(filename, ".tex")), "\n") } # ------------------------------------------------------------ # # Combine all baseline IRFs into one long dataset # ------------------------------------------------------------ # irf_main_all <- bind_rows(irfs_main) %>% mutate( shock_clean = recode( shock, ar_meu_innov = "MEU", fu_uncertainty_innov = "FU", fu_riskpremium_innov = "VRP" ), outcome_clean = recode( outcome, ip_yoy = "IP growth", inflation = "Inflation", unemployment = "Unemployment", capital_goods_growth = "Capital goods", nfc_growth = "NFC loans", ciss_std = "CISS", ciss = "CISS", credit_spreads_std = "Credit spreads", credit_spreads = "Credit spreads", d_credit_spreads = "Credit spreads" ), coef_se = fmt_coef_se(coef, se, pval), coef_star = fmt_coef(coef, pval) ) write_csv(irf_main_all, file.path(TAB_MAIN, "irf_main_all_long.csv")) openxlsx::write.xlsx( irf_main_all, file.path(TAB_MAIN, "irf_main_all_long.xlsx"), overwrite = TRUE ) # ------------------------------------------------------------ # # 1. Compact coefficient table at key horizons # ------------------------------------------------------------ # key_horizons <- c(0, 6, 12, 18, 24) irf_key_horizons <- irf_main_all %>% filter(horizon %in% key_horizons) %>% select(outcome_clean, shock_clean, horizon, coef_se) %>% mutate(horizon = paste0("h=", horizon)) %>% tidyr::pivot_wider( names_from = horizon, values_from = coef_se ) %>% arrange( factor(outcome_clean, levels = c( "IP growth", "Inflation", "Unemployment", "Capital goods", "NFC loans", "CISS", "Credit spreads" )), factor(shock_clean, levels = c("MEU", "FU", "VRP")) ) %>% rename( Outcome = outcome_clean, Shock = shock_clean ) write_csv( irf_key_horizons, file.path(TAB_MAIN, "irf_key_horizons_main.csv") ) write_tex( irf_key_horizons, filename = "irf_key_horizons_main", caption = "Baseline LP responses at selected horizons", label = "tab:irf_key_horizons_main" ) # ------------------------------------------------------------ # # 2. Peak-response summary # For each outcome-shock pair, report largest absolute response. # ------------------------------------------------------------ # irf_peak_summary <- irf_main_all %>% group_by(outcome_clean, shock_clean) %>% slice_max(order_by = abs(coef), n = 1, with_ties = FALSE) %>% ungroup() %>% transmute( Outcome = outcome_clean, Shock = shock_clean, `Peak horizon` = horizon, `Peak response` = round(coef, 3), `SE` = round(se, 3), `p-value` = round(pval, 3), `90% significant` = ifelse(sig90, "Yes", "No") ) %>% arrange( factor(Outcome, levels = c( "IP growth", "Inflation", "Unemployment", "Capital goods", "NFC loans", "CISS", "Credit spreads" )), factor(Shock, levels = c("MEU", "FU", "VRP")) ) write_csv( irf_peak_summary, file.path(TAB_MAIN, "irf_peak_summary_main.csv") ) write_tex( irf_peak_summary, filename = "irf_peak_summary_main", caption = "Peak baseline LP responses by outcome and shock", label = "tab:irf_peak_summary_main" ) # ------------------------------------------------------------ # # 3. Full IRF grids by outcome # One LaTeX file per outcome. Better than one gigantic table. # ------------------------------------------------------------ # full_grid_dir <- file.path(TEX_MAIN, "full_irf_grids") dir.create(full_grid_dir, recursive = TRUE, showWarnings = FALSE) for (yy in unique(irf_main_all$outcome_clean)) { grid_y <- irf_main_all %>% filter(outcome_clean == yy) %>% select(shock_clean, horizon, coef_star) %>% mutate(horizon = paste0("h", horizon)) %>% tidyr::pivot_wider( names_from = horizon, values_from = coef_star ) %>% arrange(factor(shock_clean, levels = c("MEU", "FU", "VRP"))) %>% rename(Shock = shock_clean) safe_name <- yy %>% stringr::str_to_lower() %>% stringr::str_replace_all("[^a-z0-9]+", "_") %>% stringr::str_replace_all("_$", "") write_csv( grid_y, file.path(TAB_MAIN, paste0("full_irf_grid_", safe_name, ".csv")) ) tex <- knitr::kable( grid_y, format = "latex", booktabs = TRUE, caption = paste0("Full baseline LP IRF grid: ", yy), label = paste0("tab:full_irf_grid_", safe_name), escape = FALSE ) writeLines(tex, file.path(full_grid_dir, paste0("full_irf_grid_", safe_name, ".tex"))) cat("Saved full grid:", file.path(full_grid_dir, paste0("full_irf_grid_", safe_name, ".tex")), "\n") } cat("\nIRF appendix tables complete.\n") cat("Compact tables saved in:", TEX_MAIN, "\n") cat("Full grids saved in:", full_grid_dir, "\n") # ============================================================ # # ANNEX: PCA ROBUSTNESS — 3x3 MACRO DASHBOARD # # ============================================================ # cat("\n--- Annex: Building PCA-CES 3x3 macro dashboard ---\n") # PCA-version shock ordering and row labels shocks_ordered_pca <- c("MEU_PCA_CES_innov", "fu_uncertainty_innov", "fu_riskpremium_innov") shock_row_labels_pca <- c( MEU_PCA_CES_innov = "MU shock (PCA-CES)", fu_uncertainty_innov = "FU shock", fu_riskpremium_innov = "VRP shock" ) # Parameterized dashboard builder (mirrors build_3x3_dashboard) build_3x3_dashboard_custom <- function(irf_list, outcomes, title_text, filename, fig_dir, shocks_ordered, shock_row_labels) { all_panels <- list() for (s in shocks_ordered) { row_panels <- build_shock_row( irf_list = irf_list, shock = s, outcomes = outcomes, labels = short_labels, ylabels = short_ylabels, color = shock_colors[s] ) all_panels <- c(all_panels, row_panels) } row_label_grobs <- purrr::map(shocks_ordered, function(s) { ggplot() + annotate("text", x = 0.5, y = 0.5, label = shock_row_labels[s], angle = 90, fontface = "bold", size = 3.2, color = shock_colors[s]) + theme_void() + theme(plot.margin = margin(0, 2, 0, 2)) }) full_layout <- list() for (i in seq_along(shocks_ordered)) { full_layout <- c( full_layout, list(row_label_grobs[[i]]), all_panels[((i-1)*length(outcomes)+1):(i*length(outcomes))] ) } n_cols <- length(outcomes) + 1 col_widths <- c(0.12, rep(1, length(outcomes))) p_final <- patchwork::wrap_plots(full_layout, ncol = n_cols, widths = col_widths) + patchwork::plot_annotation( title = title_text, subtitle = "PCA-CES robustness LP. NW-HAC. Shaded: 68% (dark) and 90% (light) CI. \u25CF = sig. at 90%.", theme = theme( plot.title = element_text(face = "bold", size = 11), plot.subtitle = element_text(size = 8.5, color = "grey35") ) ) save_fig(p_final, filename, dir = fig_dir, width = 11, height = 9) invisible(p_final) } build_3x3_dashboard_custom( irf_list = irfs_pca, outcomes = MACRO_OUTCOMES, title_text = "PCA-CES Robustness: Three Uncertainty Shocks \u00D7 Macroeconomic Outcomes", filename = "annex_dashboard_3x3_macro_pca.pdf", fig_dir = FIG_PCA, shocks_ordered = shocks_ordered_pca, shock_row_labels = shock_row_labels_pca )