# ============================================================================= # GeoRisk Data Lab 1 — Local Projections and Geopolitical Shocks # Eric Gabin Kilama · FERDI · CERDI # ============================================================================= # # References: # Jordà, Ò. (2005). "Estimation and Inference of Impulse Responses by # Local Projections." American Economic Review, 95(1), 161-182. # Caldara, D. & Iacoviello, M. (2022). "Measuring Geopolitical Risk." # American Economic Review, 112(4), 1194-1225. # # Data: GPR Index — https://www.matteoiacoviello.com/gpr.htm # ============================================================================= library(readxl) library(fixest) library(ggplot2) library(data.table) # --- 1. Load GPR data -------------------------------------------------------- # Download from https://www.matteoiacoviello.com/gpr_files/data_gpr_export.xls gpr_raw <- read_excel("data_gpr_export.xls", sheet = 1) gpr <- data.table(gpr_raw) gpr[, date := as.Date(paste0(year, "-", month, "-01"))] # Keep monthly global GPR index gpr <- gpr[, .(date, gpr = GPRD)] # --- 2. Load macro outcomes -------------------------------------------------- # Option A: FRED via fredr package # install.packages("fredr") # library(fredr) # fredr_set_key("YOUR_FRED_API_KEY") # sp500 <- fredr(series_id = "SP500", frequency = "m") # vix <- fredr(series_id = "VIXCLS", frequency = "m") # oil <- fredr(series_id = "WTISPLC", frequency = "m") # eurusd <- fredr(series_id = "DEXUSEU", frequency = "m") # Option B: CSV from FRED (download manually) sp500 <- fread("sp500_monthly.csv") vix <- fread("vix_monthly.csv") oil <- fread("wti_monthly.csv") eurusd <- fread("eurusd_monthly.csv") # Merge all series macro <- merge(gpr, sp500, by = "date", all.x = TRUE) macro <- merge(macro, vix, by = "date", all.x = TRUE) macro <- merge(macro, oil, by = "date", all.x = TRUE) macro <- merge(macro, eurusd, by = "date", all.x = TRUE) # Compute log-returns for SP500 and oil macro[, sp500_ret := c(NA, diff(log(sp500))) * 100] macro[, oil_ret := c(NA, diff(log(oil))) * 100] macro[, vix_chg := c(NA, diff(vix))] macro[, eurusd_chg := c(NA, diff(eurusd)) * 100] # Standardize GPR shock (1 SD increase) macro[, gpr_shock := (gpr - mean(gpr, na.rm = TRUE)) / sd(gpr, na.rm = TRUE)] # --- 3. Local Projections (Jordà 2005) -------------------------------------- # The key idea: for each horizon h, run a separate OLS regression: # y_{t+h} = alpha_h + beta_h * GPR_t + Gamma_h * X_t + epsilon_{t+h} # # beta_h traces out the impulse response function (IRF) H <- 12 # maximum horizon in months p <- 4 # lags to include as controls outcomes <- c("sp500_ret", "vix_chg", "oil_ret", "eurusd_chg") outcome_labels <- c("S&P 500 Return (%)", "VIX Change (pts)", "Oil Return (%)", "EUR/USD Change (%)") lp_results <- list() for (dep_var in outcomes) { irf <- numeric(H + 1) se <- numeric(H + 1) ci90_lo <- numeric(H + 1) ci90_hi <- numeric(H + 1) ci68_lo <- numeric(H + 1) ci68_hi <- numeric(H + 1) for (h in 0:H) { # Create forward variable y_{t+h} macro[, y_forward := shift(get(dep_var), n = -h, type = "lead")] # Controls: p lags of dependent variable + p lags of GPR lag_formula <- paste0( "y_forward ~ gpr_shock", paste0(" + l(", dep_var, ", ", 1:p, ")", collapse = ""), paste0(" + l(gpr_shock, ", 1:p, ")", collapse = "") ) # Estimate with Newey-West HAC standard errors (bandwidth = h+1) est <- feols(as.formula(lag_formula), data = macro, vcov = NW(h + 1)) irf[h + 1] <- coef(est)["gpr_shock"] se[h + 1] <- sqrt(vcov(est)["gpr_shock", "gpr_shock"]) ci90_lo[h + 1] <- irf[h + 1] - 1.645 * se[h + 1] ci90_hi[h + 1] <- irf[h + 1] + 1.645 * se[h + 1] ci68_lo[h + 1] <- irf[h + 1] - 1.000 * se[h + 1] ci68_hi[h + 1] <- irf[h + 1] + 1.000 * se[h + 1] } lp_results[[dep_var]] <- data.table( horizon = 0:H, irf, se, ci90_lo, ci90_hi, ci68_lo, ci68_hi ) } # --- 4. Plot IRFs ------------------------------------------------------------ plot_irf <- function(dep_var, title) { dt <- lp_results[[dep_var]] ggplot(dt, aes(x = horizon, y = irf)) + geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") + geom_ribbon(aes(ymin = ci90_lo, ymax = ci90_hi), fill = "#2C5F8A", alpha = 0.15) + geom_ribbon(aes(ymin = ci68_lo, ymax = ci68_hi), fill = "#2C5F8A", alpha = 0.30) + geom_line(color = "#1B2A4A", linewidth = 1.1) + geom_point(color = "#1B2A4A", size = 2) + labs( title = title, subtitle = "Response to 1 SD GPR shock (Jordà LP, Newey-West HAC)", x = "Horizon (months)", y = "Response" ) + scale_x_continuous(breaks = 0:H) + theme_minimal(base_family = "Helvetica Neue") + theme( plot.title = element_text(face = "bold", size = 14), plot.subtitle = element_text(color = "grey50", size = 10), panel.grid.minor = element_blank() ) } plots <- Map(plot_irf, outcomes, outcome_labels) # Arrange 2x2 panel library(patchwork) (plots[[1]] + plots[[2]]) / (plots[[3]] + plots[[4]]) + plot_annotation( title = "GeoRisk Data Lab 1 — Local Projections", subtitle = "Impulse responses of macro-financial variables to GPR shocks", caption = "Source: Caldara-Iacoviello GPR Index. Method: Jordà (2005) LP.\nEric Gabin Kilama · FERDI · CERDI", theme = theme( plot.title = element_text(face = "bold", size = 16), plot.caption = element_text(color = "grey50", hjust = 0) ) ) ggsave("lab1_irf_panel.pdf", width = 12, height = 8) ggsave("lab1_irf_panel.png", width = 12, height = 8, dpi = 300) # --- 5. State-dependent LPs (extension) ------------------------------------- # Jordà (2005) extension: IRF conditional on state of the economy # s_t = recession indicator (e.g., NBER or output gap) # macro[, recession := ifelse(output_gap < 0, 1, 0)] # # for (h in 0:H) { # macro[, y_forward := shift(get(dep_var), n = -h, type = "lead")] # est <- feols(y_forward ~ gpr_shock * recession + ..., data = macro) # } cat("\n=== Lab 1 complete ===\n") cat("IRFs saved to lab1_irf_panel.pdf\n") cat("Key finding: GPR shocks depress equity returns and increase volatility,\n") cat("with effects peaking at 2-4 months and dissipating by 8-12 months.\n")