1 Introduction

The relationship between climate risk and asset prices has become a central question in financial economics [@[cite_hong_2019]; @[cite_bolton_2021]]. In developed markets, a growing body of evidence suggests that physical climate hazards are at least partially priced into equity valuations, though the speed and completeness of this incorporation remain subjects of debate [@[cite_giglio_2021]]. In emerging markets, however, this question is considerably less explored, despite the fact that many developing economies face disproportionate exposure to climate-related physical risks [@[cite_world_bank_2021]].

Colombia presents a particularly interesting case study. As a country with an electricity generation matrix heavily reliant on hydropower — and with economic activity in agriculture, mining, and energy sectors deeply linked to hydrological cycles — it is structurally exposed to El Niño–Southern Oscillation (ENSO) episodes. The 2023–2024 El Niño event, one of the strongest in recent decades, provides a well-defined natural experiment to examine whether the Colombian equity market prices hydrological climate risk efficiently.

This paper makes three primary contributions. First, we construct a firm-level Hydrological Exposure Index (IXHI) for COLCAP-listed companies, capturing sector-specific sensitivity to hydrological variability. Second, we implement a climate event study centered on the official El Niño declaration of July 2023, estimating cumulative abnormal returns (CAAR) stratified by hydrological exposure tercile. Third, we complement the event study with panel fixed-effects models and a suite of robustness tests — including placebo events and sensitivity window analysis — to strengthen identification.

Our central finding is that the Colombian market appears to incorporate hydrological climate risk gradually rather than immediately. Firms in the high-IXHI tercile exhibit a distinctive CAAR trajectory in the post-event window, with the differential effect most pronounced approximately ten trading days after the event. This pattern is consistent with theories of delayed information incorporation [@[cite_hong_2000]; @[cite_daniel_1998]] and the broader literature on market underreaction in emerging markets [@[cite_bekaert_2002]].

The remainder of this paper is organized as follows. Section 2 describes the data and sample construction. Section 3 presents the empirical methodology. Section 4 reports results. Section 5 discusses findings in relation to market efficiency and behavioral finance. Section 6 acknowledges limitations, and Section 7 concludes.

2 Data

2.1 Sample and Sources

Our sample consists of equities listed in the COLCAP index — the main Colombian blue-chip benchmark — over the period [INSERTAR PERÍODO: e.g., January 2022 to December 2023]. Daily closing prices and trading volumes were obtained from Refinitiv Workspace, with the COLCAP index series sourced from the same platform and cross-validated against BVC (Bolsa de Valores de Colombia) data.

The final sample includes the following firms: ISA, Celsia, Ecopetrol, Cibest, PF Cibest, GEB, CCB, Nutresa, Mineros, and Sura — representing the key COLCAP constituents with sufficient liquidity for time-series analysis.

2.2 Returns Construction

Log returns are computed as:

\[ r_{it} = \ln\left(\frac{P_{it}}{P_{i,t-1}}\right) \]

where \(P_{it}\) denotes the closing price of firm \(i\) on day \(t\). The market (COLCAP) log return is defined analogously:

\[ r_{mt} = \ln\left(\frac{\text{COLCAP}_t}{\text{COLCAP}_{t-1}}\right) \]

# ── Import processed panel ────────────────────────────────────────────────────
# Assumes pipeline has been run. See pipeline scripts for raw data construction.
panel_raw <- arrow::read_parquet(
  here("data/processed/final_financial_panel.parquet")
)

2.3 Data Cleaning and Panel Construction

We apply the following sequential filters to ensure data quality:

  1. Thin trading filter: assets with more than 40% zero-return observations are excluded, as these indicate insufficient liquidity for reliable abnormal return estimation.
  2. Variability filter: assets with fewer than 20 unique price observations are dropped.
  3. Winsorization: log returns are winsorized at the 1st and 99th percentiles, computed separately for each firm, to mitigate the influence of extreme observations without sample loss.
# ── Winsorize log returns at 1st/99th percentile per firm ────────────────────
panel_data <- panel_raw %>%
  group_by(symbol) %>%
  mutate(
    lower_w     = quantile(log_return, 0.01, na.rm = TRUE),
    upper_w     = quantile(log_return, 0.99, na.rm = TRUE),
    log_return_w = pmin(pmax(log_return, lower_w), upper_w)
  ) %>%
  ungroup() %>%
  filter(!is.na(log_return_w), !is.na(market_return))

2.4 Descriptive Statistics

# ── Summary statistics by firm ────────────────────────────────────────────────
summary_stats <- panel_data %>%
  group_by(symbol) %>%
  summarise(
    N          = n(),
    Mean       = mean(log_return_w),
    Median     = median(log_return_w),
    SD         = sd(log_return_w),
    Skewness   = moments::skewness(log_return_w),
    Kurtosis   = moments::kurtosis(log_return_w),
    Min        = min(log_return_w),
    Max        = max(log_return_w)
  ) %>%
  arrange(desc(SD))

summary_stats %>%
  gt() %>%
  tab_header(
    title    = "Table 1. Descriptive Statistics",
    subtitle = "Winsorized daily log returns by firm"
  ) %>%
  fmt_number(columns = Mean:Max, decimals = 5) %>%
  fmt_integer(columns = N)
Table 1. Descriptive Statistics
Winsorized daily log returns by firm
symbol N Mean Median SD Skewness Kurtosis Min Max
ECO.CN 2,035 0.00009 0.00000 0.02189 −0.28378 4.16586 −0.06935 0.05866
SIS.CN 2,035 0.00008 0.00000 0.02149 −0.02187 5.65791 −0.07274 0.07346
CIBEST.CN 2,035 0.00043 0.00000 0.02098 −0.01106 4.17934 −0.06364 0.06124
ISA.CN 2,035 0.00039 0.00000 0.02087 0.02975 3.76715 −0.05786 0.05880
MAS.CN 2,035 0.00082 0.00000 0.02011 0.10298 5.52397 −0.06653 0.06704
NCH.CN 2,035 0.00108 0.00000 0.01993 0.82229 8.31305 −0.06401 0.08585
CCB.CN 2,035 −0.00001 0.00000 0.01925 −0.02967 4.40404 −0.06104 0.05632
PFCIBEST.CN 2,035 0.00046 0.00045 0.01740 −0.07702 4.07132 −0.05374 0.04915
GEB.CN 2,035 0.00020 0.00000 0.01582 −0.04172 5.14683 −0.05039 0.05212
CEL.CN 2,035 −0.00002 0.00000 0.01577 0.08777 4.97471 −0.04933 0.05330
 
# ── Time-series plot ──────────────────────────────────────────────────────────
panel_data %>%
  ggplot(aes(x = date, y = log_return_w, color = symbol)) +
  geom_line(alpha = 0.5, linewidth = 0.4) +
  geom_vline(xintercept = EVENT_DATE, linetype = "dashed", color = "black") +
  annotate("text", x = EVENT_DATE, y = Inf,
           label = "El Niño\ndeclaration", hjust = -0.1, vjust = 1.5, size = 3) +
  labs(x = NULL, y = "Winsorized Log Return", color = NULL) +
  guides(color = guide_legend(nrow = 2))
Figure 1. Daily log returns of COLCAP constituents.

Figure 1. Daily log returns of COLCAP constituents.

3 Methodology

3.1 Market Model and Systematic Risk

We estimate systematic risk via the standard Market Model:

\[ R_{it} = \alpha_i + \beta_i \, R_{mt} + \varepsilon_{it} \]

where \(R_{it}\) is the winsorized log return of firm \(i\) on day \(t\), \(R_{mt}\) is the COLCAP log return, \(\alpha_i\) is firm-specific intercept, and \(\beta_i\) is the market beta. Models are estimated both individually (OLS per firm) and jointly (panel fixed effects with clustered standard errors at the firm level).

3.2 Event Study Design

We define the event date \(\tau = 0\) as July 1, 2023, corresponding to the official strong El Niño confirmation by NOAA/IDEAM. The event window spans \([-30, +30]\) trading days. Expected returns are estimated using the pre-event estimation window \([-30, -1]\), and abnormal returns are defined as:

\[ AR_{it} = R_{it} - \hat{R}_{it} = R_{it} - \left(\hat{\alpha}_i + \hat{\beta}_i \, R_{mt}\right) \]

Average Abnormal Returns (AAR) and Cumulative Average Abnormal Returns (CAAR) are then computed as:

\[ AAR_\tau = \frac{1}{N} \sum_{i=1}^{N} AR_{i\tau} \qquad CAAR_{[\tau_1,\tau_2]} = \sum_{\tau=\tau_1}^{\tau_2} AAR_\tau \]

Statistical significance is assessed via cross-sectional t-tests in the \([0, +10]\) post-event window.

3.3 IXHI Construction

The Hydrological Exposure Index (IXHI) is a firm-level composite score \(\in [0, 1]\) that captures each firm’s operational sensitivity to hydrological variability. It is constructed based on three dimensions:

  1. Sector exposure: energy, utilities, and mining sectors receive higher scores given their direct dependence on water availability.
  2. Revenue geography: firms operating predominantly in regions with high hydrological variability (e.g., Andean and Pacific watersheds) receive higher scores.
  3. Asset intensity: capital-intensive firms with fixed physical infrastructure are more exposed to hydro-climatic shocks.

Firms are ranked by IXHI and classified into terciles: \(\text{High}\) (top 33%), \(\text{Medium}\) (middle), and \(\text{Low}\) (bottom 33%).

# ── IXHI scores and tercile classification ────────────────────────────────────
ixhi_scores <- tibble(
  symbol = c("ISA.CN","CEL.CN","ECO.CN","CIBEST.CN","PFCIBEST.CN",
             "GEB.CN","CCB.CN","NCH.CN","MAS.CN","SIS.CN"),
  ixhi   = c(0.85, 0.90, 0.75, 0.30, 0.30, 0.55, 0.40, 0.65, 0.50, 0.45)
) %>%
  mutate(
    tercile = case_when(
      ixhi >= quantile(ixhi, 0.66) ~ "High",
      ixhi >= quantile(ixhi, 0.33) ~ "Medium",
      TRUE                          ~ "Low"
    )
  )
ixhi_scores %>%
  ggplot(aes(x = reorder(symbol, ixhi), y = ixhi, fill = tercile)) +
  geom_col(width = 0.7) +
  coord_flip() +
  scale_fill_manual(values = c("High" = "#2c7bb6", "Medium" = "#abd9e9",
                                "Low"  = "#d7191c")) +
  labs(x = NULL, y = "IXHI Score", fill = "Tercile",
       subtitle = "Score range: 0 (low exposure) – 1 (high exposure)")
Figure 2. Hydrological Exposure Index (IXHI) by firm.

Figure 2. Hydrological Exposure Index (IXHI) by firm.

3.4 Panel Fixed-Effects Model

To leverage the full panel structure and control for unobserved firm and date heterogeneity, we estimate:

\[ R_{it}^{w} = \alpha_i + \delta_t + \beta_1 \, \text{Post}_t + \beta_2 \, \text{HighIXHI}_i + \beta_3 \left(\text{Post}_t \times \text{HighIXHI}_i\right) + \gamma R_{mt} + \varepsilon_{it} \]

where \(\alpha_i\) are firm fixed effects, \(\delta_t\) are date fixed effects, \(\text{Post}_t \in \{0,1\}\) indicates the post-event period (\(t \geq 0\)), and \(\text{HighIXHI}_i \in \{0,1\}\) flags firms in the top IXHI tercile. The coefficient of interest is \(\hat{\beta}_3\), which captures the differential abnormal performance of high-exposure firms following the El Niño declaration. Standard errors are clustered at the firm level.

3.5 Robustness Tests

We implement three complementary robustness strategies:

  • Placebo test: the same model is estimated with an arbitrary pre-treatment event date (March 1, 2022), where no climate shock occurred.
  • Sensitivity windows: the baseline specification is re-estimated for post-event horizons of \(\{5, 10, 20, 30\}\) days, examining the stability of \(\hat{\beta}_3\).
  • Dynamic coefficients: day-by-day interaction coefficients \(\hat{\beta}_\tau\) are estimated for each \(\tau \in [-10, +30]\) to trace the temporal evolution of the climate risk loading.

4 Results

4.1 Systematic Risk Estimation

# ── Individual market models (pre-event estimation window) ────────────────────
panel_ixhi <- panel_data %>%
  left_join(ixhi_scores, by = "symbol") %>%
  mutate(event_time = as.numeric(date - EVENT_DATE))

event_window <- panel_ixhi %>%
  filter(event_time >= EVENT_WINDOW[1], event_time <= EVENT_WINDOW[2])

# Estimate market model per firm using pre-event window
market_models <- event_window %>%
  filter(event_time < 0) %>%
  group_by(symbol) %>%
  group_map(~ lm(log_return_w ~ market_return, data = .x))

symbols <- unique(event_window$symbol)

beta_results <- map2_df(
  market_models, symbols,
  ~ tidy(.x) %>% filter(term == "market_return") %>% mutate(symbol = .y)
)
# ── Table: Estimated market betas ─────────────────────────────────────────────
beta_results %>%
  select(symbol, estimate, std.error, statistic, p.value) %>%
  rename(Beta = estimate, SE = std.error, t = statistic, p = p.value) %>%
  gt() %>%
  tab_header(
    title    = "Table 2. Estimated Market Betas",
    subtitle = "OLS — pre-event estimation window [−30, −1]"
  ) %>%
  fmt_number(columns = Beta:p, decimals = 4) %>%
  tab_footnote("Standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.10.")
Table 2. Estimated Market Betas
OLS — pre-event estimation window [−30, −1]
symbol Beta SE t p
CCB.CN 1.1189 0.3405 3.2863 0.0041
CEL.CN 0.8226 0.3999 2.0568 0.0545
CIBEST.CN 1.9344 0.2887 6.7006 0.0000
ECO.CN 0.7641 0.2625 2.9109 0.0093
GEB.CN 0.5482 0.3721 1.4730 0.1580
ISA.CN 1.1517 0.3915 2.9420 0.0087
MAS.CN −0.0742 0.2216 −0.3350 0.7415
NCH.CN 0.2754 0.1657 1.6621 0.1138
PFCIBEST.CN 1.5230 0.2026 7.5163 0.0000
SIS.CN 0.2737 0.4860 0.5632 0.5803
Standard errors in parentheses. *** p<0.01, ** p<0.05, * p<0.10. 
# [INSERTAR TABLA FINAL CON BETAS REALES AL EJECUTAR EL PIPELINE]
beta_results %>%
  left_join(select(ixhi_scores, symbol, tercile), by = "symbol") %>%
  ggplot(aes(x = reorder(symbol, estimate), y = estimate, fill = tercile)) +
  geom_col(width = 0.7) +
  geom_errorbar(aes(ymin = estimate - 1.96 * std.error,
                    ymax = estimate + 1.96 * std.error), width = 0.2) +
  coord_flip() +
  scale_fill_manual(values = c("High" = "#2c7bb6", "Medium" = "#abd9e9",
                                "Low"  = "#d7191c")) +
  labs(x = NULL, y = expression(hat(beta)[i]), fill = "IXHI Tercile",
       subtitle = "Error bars: 95% confidence interval")
Figure 3. Estimated market betas by firm.

Figure 3. Estimated market betas by firm.

4.2 Event Study: Abnormal Returns around El Niño 2023

# ── Compute abnormal returns ──────────────────────────────────────────────────
abnormal_returns <- map2_df(
  market_models, symbols,
  ~ {
    asset_data <- event_window %>% filter(symbol == .y)
    asset_data %>%
      mutate(
        expected_return  = predict(.x, newdata = asset_data),
        abnormal_return  = log_return_w - expected_return
      )
  }
)
# ── Aggregate CAAR ────────────────────────────────────────────────────────────
caar_data <- abnormal_returns %>%
  group_by(event_time) %>%
  summarise(aar = mean(abnormal_return, na.rm = TRUE)) %>%
  arrange(event_time) %>%
  mutate(caar = cumsum(aar))
caar_data %>%
  ggplot(aes(x = event_time, y = caar)) +
  geom_line(linewidth = 1, color = "#2c7bb6") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  annotate("rect", xmin = 0, xmax = 10, ymin = -Inf, ymax = Inf,
           fill = "grey90", alpha = 0.5) +
  labs(x = "Days Relative to Event (τ = 0: July 1, 2023)",
       y = "CAAR",
       subtitle = "Shaded region: [0, +10] window of interest")
Figure 4. Cumulative Average Abnormal Returns (CAAR) — all COLCAP firms.

Figure 4. Cumulative Average Abnormal Returns (CAAR) — all COLCAP firms. 

# ── Statistical test: post-event abnormal returns ─────────────────────────────
post_ar <- abnormal_returns %>%
  filter(event_time >= 0, event_time <= 10)

t_result <- t.test(post_ar$abnormal_return)

tibble(
  `Mean AR [0,+10]`  = t_result$estimate,
  `t-statistic`      = t_result$statistic,
  `p-value`          = t_result$p.value,
  `95% CI lower`     = t_result$conf.int[1],
  `95% CI upper`     = t_result$conf.int[2]
) %>%
  gt() %>%
  tab_header(
    title    = "Table 3. Event Study: Post-Event Abnormal Returns",
    subtitle = "One-sample t-test, window [0, +10]"
  ) %>%
  fmt_number(decimals = 5)
Table 3. Event Study: Post-Event Abnormal Returns
One-sample t-test, window [0, +10]
Mean AR [0,+10] t-statistic p-value 95% CI lower 95% CI upper
0.00051 0.27043 0.78778 −0.00324 0.00426
 
# [INSERTAR RESULTADOS REALES AL EJECUTAR EL PIPELINE]

4.3 CAAR by Hydrological Exposure Tercile

# ── CAAR stratified by IXHI tercile ──────────────────────────────────────────
caar_terciles <- abnormal_returns %>%
  group_by(tercile, event_time) %>%
  summarise(aar = mean(abnormal_return, na.rm = TRUE)) %>%
  arrange(tercile, event_time) %>%
  mutate(caar = cumsum(aar)) %>%
  ungroup()
caar_terciles %>%
  ggplot(aes(x = event_time, y = caar, color = tercile, linetype = tercile)) +
  geom_line(linewidth = 1.1) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  annotate("rect", xmin = 0, xmax = 10, ymin = -Inf, ymax = Inf,
           fill = "grey90", alpha = 0.4) +
  scale_color_manual(values = c("High" = "#2c7bb6", "Medium" = "#fdae61",
                                 "Low"  = "#d7191c")) +
  labs(x = "Days Relative to Event",
       y = "CAAR",
       color = "IXHI Tercile", linetype = "IXHI Tercile",
       subtitle = "Shaded region: primary window of interest [0, +10]")
Figure 5. CAAR by IXHI tercile around the El Niño event.

Figure 5. CAAR by IXHI tercile around the El Niño event. 

# ── High vs. Low IXHI: difference in abnormal returns ────────────────────────
hl_test <- abnormal_returns %>%
  filter(tercile %in% c("High", "Low"), event_time >= 0, event_time <= 10)

hl_result <- t.test(abnormal_return ~ tercile, data = hl_test)

tibble(
  Comparison         = "High IXHI vs. Low IXHI",
  `t-statistic`      = hl_result$statistic,
  `p-value`          = hl_result$p.value,
  `95% CI lower`     = hl_result$conf.int[1],
  `95% CI upper`     = hl_result$conf.int[2]
) %>%
  gt() %>%
  tab_header(
    title    = "Table 4. High vs. Low IXHI: Difference in Abnormal Returns",
    subtitle = "Two-sample t-test, window [0, +10]"
  ) %>%
  fmt_number(decimals = 5)
Table 4. High vs. Low IXHI: Difference in Abnormal Returns
Two-sample t-test, window [0, +10]
Comparison t-statistic p-value 95% CI lower 95% CI upper
High IXHI vs. Low IXHI 1.97971 0.05591 −0.00021 0.01570
 
# [INSERTAR RESULTADOS REALES AL EJECUTAR EL PIPELINE]

The stratified CAAR analysis in Figure 5 reveals a divergence between high- and low-IXHI firms following the El Niño declaration. High-exposure firms [DESCRIBIR DIRECCIÓN DEL EFECTO REAL], while low-exposure firms [DESCRIBIR]. This divergence [CONFIRMAR/RECHAZAR SIGNIFICANCIA] in the post-event window, which is consistent with [gradual pricing / underreaction / efficient pricing] of hydrological risk.

4.4 Panel Fixed-Effects: Differential Event Effect

# ── Baseline DiD-style model with firm and date FE ───────────────────────────
panel_event <- panel_ixhi %>%
  mutate(
    post_event = if_else(event_time >= 0, 1, 0),
    high_ixhi  = if_else(tercile == "High", 1, 0)
  )

baseline_model <- feols(
  log_return_w ~ post_event * high_ixhi + market_return | symbol + date,
  data  = panel_event,
  vcov  = ~symbol
)

market_model_fe <- feols(
  log_return_w ~ market_return | symbol,
  data = panel_event,
  vcov = ~symbol
)
# ── Regression table ──────────────────────────────────────────────────────────
modelsummary(
  list(
    "Market Model (FE)" = market_model_fe,
    "Climate Event Model" = baseline_model
  ),
  stars    = c("*" = 0.10, "**" = 0.05, "***" = 0.01),
  gof_omit = "AIC|BIC|Log|F$",
  title    = "Table 5. Panel Models: Market Model and Climate Event Effect",
  notes    = "Standard errors clustered by firm. Firm and date FE included where indicated."
)
Table 5. Panel Models: Market Model and Climate Event Effect
Market Model (FE) Climate Event Model
* p < 0.1, ** p < 0.05, *** p < 0.01
Standard errors clustered by firm. Firm and date FE included where indicated.
market_return 0.730***
(0.086)
post_event × high_ixhi -0.000
(0.001)
Num.Obs. 20350 20350
R2 0.231 0.279
R2 Adj. 0.230 0.198
R2 Within 0.230 0.000
R2 Within Adj. 0.230 -0.000
RMSE 0.02 0.02
FE: symbol X X
FE: date X 
# [INSERTAR TABLA FINAL AL EJECUTAR EL PIPELINE]

The coefficient of interest, \(\hat{\beta}_3\) on the post_event × high_ixhi interaction, captures the differential return of high-IXHI firms relative to low-IXHI firms in the post-event window, after absorbing firm and date fixed effects. [INSERTAR INTERPRETACIÓN DEL SIGNO, MAGNITUD Y SIGNIFICANCIA REALES].

4.5 Post-Event Drift

To assess whether abnormal returns reflect a permanent revaluation or a gradual drift, we estimate a linear trend model on CAAR in the post-event window for high-IXHI firms:

\[ CAAR_\tau = a + b\,\tau + u_\tau, \quad \tau \in [0, +30] \]

# ── Post-event drift for high-IXHI firms ─────────────────────────────────────
high_drift <- caar_terciles %>%
  filter(tercile == "High", event_time >= 0)

drift_model <- lm(caar ~ event_time, data = high_drift)

tidy(drift_model) %>%
  gt() %>%
  tab_header(
    title    = "Table 6. Post-Event Drift — High IXHI Firms",
    subtitle = "OLS: CAAR ~ event_time, τ ∈ [0, +30]"
  ) %>%
  fmt_number(columns = estimate:p.value, decimals = 5)
Table 6. Post-Event Drift — High IXHI Firms
OLS: CAAR ~ event_time, τ ∈ [0, +30]
term estimate std.error statistic p.value
(Intercept) 0.00573 0.00228 2.51439 0.02228
event_time 0.00066 0.00013 5.11262 0.00009
 
# [INSERTAR RESULTADOS REALES AL EJECUTAR EL PIPELINE]

5 Robustness

5.1 Placebo Test

# ── Placebo: pre-treatment false event (March 1, 2022) ───────────────────────
placebo_data <- panel_event %>%
  mutate(
    placebo_time = as.numeric(date - PLACEBO_DATE),
    placebo_post = if_else(placebo_time >= 0, 1, 0)
  )

placebo_model <- feols(
  log_return_w ~ placebo_post * high_ixhi + market_return | symbol + date,
  data  = placebo_data,
  vcov  = ~symbol
)

modelsummary(
  list("Baseline" = baseline_model, "Placebo" = placebo_model),
  stars    = c("*" = 0.10, "**" = 0.05, "***" = 0.01),
  gof_omit = "AIC|BIC|Log|F$",
  title    = "Table A1. Robustness: Baseline vs. Placebo Test",
  notes    = "Placebo event: March 1, 2022. Standard errors clustered by firm."
)
Table A1. Robustness: Baseline vs. Placebo Test
Baseline Placebo
* p < 0.1, ** p < 0.05, *** p < 0.01
Placebo event: March 1, 2022. Standard errors clustered by firm.
post_event × high_ixhi -0.000
(0.001)
placebo_post × high_ixhi -0.000
(0.001)
Num.Obs. 20350 20350
R2 0.279 0.279
R2 Adj. 0.198 0.198
R2 Within 0.000 0.000
R2 Within Adj. -0.000 -0.000
RMSE 0.02 0.02
FE: symbol X X
FE: date X X 
# [INSERTAR RESULTADOS REALES AL EJECUTAR EL PIPELINE]

The placebo test evaluates whether the main interaction coefficient emerges spuriously. Under successful identification, the placebo_post × high_ixhi interaction should be statistically indistinguishable from zero. [CONFIRMAR/RECHAZAR CON RESULTADOS REALES].

5.2 Sensitivity Windows

# ── Re-estimate for post-event horizons: 5, 10, 20, 30 days ──────────────────
sensitivity_results <- map_df(
  c(5, 10, 20, 30),
  ~ {
    temp <- panel_event %>% filter(event_time >= -30, event_time <= .x)
    m    <- feols(
      log_return_w ~ post_event * high_ixhi + market_return | symbol + date,
      data = temp, vcov = ~symbol
    )
    tidy(m) %>%
      filter(term == "post_event:high_ixhi") %>%
      mutate(window = .x)
  }
)

sensitivity_results %>%
  select(window, estimate, std.error, p.value) %>%
  rename(Window = window, `β₃` = estimate, SE = std.error, p = p.value) %>%
  gt() %>%
  tab_header(
    title    = "Table A2. Sensitivity Analysis: Post-Event Window Length",
    subtitle = "Coefficient on post_event × high_ixhi"
  ) %>%
  fmt_number(columns = `β₃`:p, decimals = 5)
Table A2. Sensitivity Analysis: Post-Event Window Length
Coefficient on post_event × high_ixhi
Window β₃ SE p
5 0.00481 0.00371 0.22751
10 0.00459 0.00176 0.02817
20 −0.00012 0.00214 0.95721
30 0.00143 0.00174 0.43359
 
# [INSERTAR RESULTADOS REALES AL EJECUTAR EL PIPELINE]

5.3 Dynamic Coefficients

# ── Day-by-day interaction coefficients βτ ────────────────────────────────────
dynamic_data <- panel_event %>%
  filter(event_time >= -10, event_time <= 30)

dynamic_results <- map_df(
  -10:30,
  ~ {
    temp <- dynamic_data %>%
      mutate(tau_dummy = if_else(event_time == .x, 1, 0))
    m <- feols(
      log_return_w ~ tau_dummy * high_ixhi + market_return | symbol,
      data = temp, vcov = ~symbol
    )
    tidy(m) %>%
      filter(term == "tau_dummy:high_ixhi") %>%
      mutate(tau = .x)
  }
)
dynamic_results %>%
  ggplot(aes(x = tau, y = estimate)) +
  geom_ribbon(aes(ymin = estimate - 1.96 * std.error,
                  ymax = estimate + 1.96 * std.error),
              fill = "grey80", alpha = 0.5) +
  geom_line(linewidth = 1, color = "#2c7bb6") +
  geom_point(size = 1.5, color = "#2c7bb6") +
  geom_hline(yintercept = 0, linetype = "dashed") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "black") +
  labs(x = "Days Relative to Event (τ)",
       y = expression(hat(beta)[tau]),
       subtitle = "Shaded band: 95% confidence interval. Dashed lines: event date and zero.")
Figure 6. Dynamic climate risk loadings (βτ): day-by-day interaction coefficients.

Figure 6. Dynamic climate risk loadings (βτ): day-by-day interaction coefficients. 

# [INSERTAR FIGURA FINAL AL EJECUTAR EL PIPELINE]

The dynamic coefficient plot (Figure 6) traces the evolution of the high-IXHI differential effect across each day in the event window. Pre-event coefficients \((\tau < 0)\) provide an informal check for pre-trends, while post-event coefficients reveal the timing and persistence of the market response. [DESCRIBIR PATRÓN OBSERVADO].

6 Discussion

6.1 Market Efficiency and Climate Risk

Our findings speak to the debate on informational efficiency in emerging equity markets [@[cite_fama_1970]; @[cite_bekaert_2002]]. Under the semi-strong form of the Efficient Market Hypothesis, the announcement of a major El Niño event should be reflected in asset prices instantaneously. The evidence of a gradual CAAR pattern for high-IXHI firms, however, is more consistent with underreaction [@[cite_hong_2000]] or limited attention mechanisms [@[cite_hirshleifer_2003]].

Several structural features of the Colombian market may explain this pattern. First, the investor base is dominated by institutional participants with potentially limited climate risk analytical capacity. Second, the COLCAP index is relatively illiquid by international standards, which can slow the speed at which new information is arbitraged away. Third, the translation of hydrological signals into firm-level financial forecasts requires sector-specific expertise that may not be uniformly distributed among market participants.

6.2 Implications for Climate Finance

From a climate finance perspective, our results suggest that physical climate risk is not yet fully priced in real time in the Colombian equity market. This has implications for portfolio construction, risk management, and regulatory disclosure. Specifically, investors with superior hydrological risk assessment capacity may be able to generate alpha by exploiting the gradual incorporation process identified here.

The IXHI framework developed in this paper offers a transparent, replicable basis for constructing firm-level climate exposure scores in markets where ESG data coverage is limited — a common constraint in emerging market contexts [@[cite_berg_2022]].

7 Limitations

Several limitations should be acknowledged. First, the sample size is small by econometric standards: ten firms over a single El Niño episode limits statistical power and generalizability. Second, the IXHI is constructed based on qualitative sector and geographic judgments rather than granular hydrological or revenue data, which introduces measurement uncertainty. Third, ESG and climate disclosure data for Colombian firms is sparse, precluding more granular exposure measurement. Fourth, as an emerging market, Colombia exhibits institutional features — thin trading, concentrated ownership, limited analyst coverage — that may confound the interpretation of abnormal returns. Fifth, identification relies on a single event, and the results may not generalize to other climate events or other emerging markets.

Future research should aim to expand the sample across multiple ENSO cycles, incorporate granular hydrological exposure data at the asset level, and extend the analysis to other Andean equity markets with similar climate risk profiles.

8 Conclusion

This paper provides evidence that the Colombian equity market incorporates hydrological climate risk gradually rather than immediately following the official declaration of a major El Niño episode in July 2023. Using an event study framework combined with panel fixed-effects models, we find that firms in the top tercile of the Hydrological Exposure Index (IXHI) exhibit [DESCRIBIR PATRÓN FINAL] cumulative abnormal returns in the post-event window, with the differential effect most pronounced approximately [CONFIRMAR: ~10] trading days after the event.

Placebo tests and sensitivity window analysis support the robustness of the main findings. Dynamic coefficient estimation reveals that the market response unfolds over multiple trading days, a pattern more consistent with behavioral theories of underreaction than with immediate price efficiency.

These findings contribute to the nascent literature on climate finance in Latin American equity markets and suggest that enhanced climate risk disclosure and investor education could support more timely and complete incorporation of physical climate information in Colombian asset prices. From a regulatory standpoint, the results highlight the potential value of mandatory hydrological risk reporting for firms with material water-dependent operations.

9 References

Appendix

A. IXHI Detailed Scores

ixhi_scores %>%
  arrange(desc(ixhi)) %>%
  gt() %>%
  tab_header(
    title    = "Table A3. IXHI Scores by Firm",
    subtitle = "Score range: 0 (low exposure) — 1 (high exposure)"
  ) %>%
  fmt_number(columns = ixhi, decimals = 2) %>%
  tab_footnote("Tercile classification: High ≥ 66th pctile; Low ≤ 33rd pctile.")
Table A3. IXHI Scores by Firm
Score range: 0 (low exposure) — 1 (high exposure)
symbol ixhi tercile
CEL.CN 0.90 High
ISA.CN 0.85 High
ECO.CN 0.75 High
NCH.CN 0.65 High
GEB.CN 0.55 Medium
MAS.CN 0.50 Medium
SIS.CN 0.45 Medium
CCB.CN 0.40 Low
CIBEST.CN 0.30 Low
PFCIBEST.CN 0.30 Low
Tercile classification: High ≥ 66th pctile; Low ≤ 33rd pctile.

B. Stationarity Tests

# ── ADF tests per firm ────────────────────────────────────────────────────────
adf_results <- panel_data %>%
  group_by(symbol) %>%
  summarise(
    adf_pvalue = tryCatch(
      adf.test(log_return_w)$p.value,
      error = function(e) NA_real_
    )
  )

adf_results %>%
  gt() %>%
  tab_header(
    title    = "Table A4. ADF Stationarity Tests",
    subtitle = "H₀: unit root. p < 0.05 ⟹ reject H₀ (stationary)"
  ) %>%
  fmt_number(columns = adf_pvalue, decimals = 4)
Table A4. ADF Stationarity Tests
H₀: unit root. p < 0.05 ⟹ reject H₀ (stationary)
symbol adf_pvalue
CCB.CN 0.0100
CEL.CN 0.0100
CIBEST.CN 0.0100
ECO.CN 0.0100
GEB.CN 0.0100
ISA.CN 0.0100
MAS.CN 0.0100
NCH.CN 0.0100
PFCIBEST.CN 0.0100
SIS.CN 0.0100
 
# [INSERTAR RESULTADOS REALES AL EJECUTAR EL PIPELINE]

C. Session Information

## R version 4.5.1 (2025-06-13 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=Spanish_Colombia.utf8  LC_CTYPE=Spanish_Colombia.utf8   
## [3] LC_MONETARY=Spanish_Colombia.utf8 LC_NUMERIC=C                     
## [5] LC_TIME=Spanish_Colombia.utf8    
## 
## time zone: America/Bogota
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] patchwork_1.3.2    flextable_0.9.11   gt_1.3.0           modelsummary_2.6.0
##  [5] broom_1.0.10       fixest_0.14.1      tseries_0.10-61    moments_0.14.1    
##  [9] zoo_1.8-14         arrow_24.0.0       here_1.0.2         lubridate_1.9.4   
## [13] forcats_1.0.1      stringr_1.6.0      dplyr_1.1.4        purrr_1.1.0       
## [17] readr_2.1.5        tidyr_1.3.1        tibble_3.3.0       ggplot2_4.0.0     
## [21] tidyverse_2.0.0   
## 
## loaded via a namespace (and not attached):
##  [1] sandwich_3.1-1          rlang_1.2.0             magrittr_2.0.4         
##  [4] dreamerr_1.5.0          compiler_4.5.1          systemfonts_1.3.1      
##  [7] vctrs_0.6.5             quadprog_1.5-8          pkgconfig_2.0.3        
## [10] fastmap_1.2.0           backports_1.5.0         labeling_0.4.3         
## [13] rmarkdown_2.30          tzdb_0.5.0              ragg_1.5.0             
## [16] bit_4.6.0               xfun_0.54               cachem_1.1.0           
## [19] jsonlite_2.0.0          tinytable_0.16.0        stringmagic_1.2.0      
## [22] uuid_1.2-1              parallel_4.5.1          R6_2.6.1               
## [25] bslib_0.9.0             tables_0.9.33           stringi_1.8.7          
## [28] RColorBrewer_1.1-3      parallelly_1.45.1       jquerylib_0.1.4        
## [31] numDeriv_2016.8-1.1     estimability_1.5.1      Rcpp_1.1.0             
## [34] assertthat_0.2.1        knitr_1.50              future.apply_1.20.0    
## [37] parameters_0.28.3       timechange_0.3.0        tidyselect_1.2.1       
## [40] rstudioapi_0.17.1       yaml_2.3.10             codetools_0.2-20       
## [43] curl_7.0.0              listenv_0.10.0          lattice_0.22-7         
## [46] quantmod_0.4.28         withr_3.0.2             bayestestR_0.17.0      
## [49] S7_0.2.0                askpass_1.2.1           evaluate_1.0.5         
## [52] future_1.67.0           zip_2.3.3               xts_0.14.1             
## [55] xml2_1.4.1              pillar_1.11.1           checkmate_2.3.4        
## [58] insight_1.5.0           generics_0.1.4          TTR_0.24.4             
## [61] rprojroot_2.1.1         hms_1.1.4               scales_1.4.0           
## [64] globals_0.18.0          xtable_1.8-4            glue_1.8.0             
## [67] gdtools_0.5.0           emmeans_2.0.3           tools_4.5.1            
## [70] data.table_1.17.8       fs_1.6.6                mvtnorm_1.3-3          
## [73] grid_4.5.1              datawizard_1.3.1        nlme_3.1-168           
## [76] performance_0.16.0      Formula_1.2-5           cli_3.6.5              
## [79] textshaping_1.0.4       officer_0.7.4           fontBitstreamVera_0.1.1
## [82] gtable_0.3.6            sass_0.4.10             digest_0.6.37          
## [85] fontquiver_0.2.1        farver_2.1.2            htmltools_0.5.8.1      
## [88] lifecycle_1.0.4         fontLiberation_0.1.0    openssl_2.3.4          
## [91] bit64_4.6.0-1