# ==============================================================================
# THESIS — Liquidity and the Cross-Section of Stock Returns in Germany and the effect of different market states
# SCRIPT FINAL PROPRE ET ORGANISÉ
# BLOC 1/3 — SETUP + BASELINE + TABLES 1 TO 17
# ==============================================================================

# ==============================================================================
# 0. PACKAGES
# ==============================================================================
required_packages <- c(
  "haven", "readr", "dplyr", "tidyr", "slider", "purrr",
  "broom", "lmtest", "sandwich", "writexl", "tibble"
)

installed <- rownames(installed.packages())
to_install <- setdiff(required_packages, installed)
if (length(to_install) > 0) install.packages(to_install)

invisible(lapply(required_packages, library, character.only = TRUE))

# ==============================================================================
# 1. PATHS
# ==============================================================================
PATH_MONTHLY <- "C:/Users/najwa/Downloads/germany_monthly_long_total_stock_after_screens_all.dta"
PATH_DAILY   <- "C:/Users/najwa/Downloads/germany_daily_long_total_stock_after_screens_all.dta"
PATH_SAFE    <- "C:/Users/najwa/Downloads/SAFE_Manager_Sentiment_Index_Data_Webpage_202504.csv"
OUT_DIR      <- "C:/Users/najwa/Downloads"

PATH_BASELINE_RDS          <- file.path(OUT_DIR, "baseline_sample_final.rds")
PATH_BASELINE_SENT_RDS     <- file.path(OUT_DIR, "baseline_sample_with_sentiment.rds")
PATH_BASELINE_EDGE_RDS     <- file.path(OUT_DIR, "baseline_sample_with_edge.rds")
PATH_BASELINE_EDGE_ONLYRDS <- file.path(OUT_DIR, "baseline_sample_edge_only.rds")

# ==============================================================================
# 2. HELPER FUNCTIONS
# ==============================================================================
round_numeric <- function(df, digits = 6) {
  df %>% mutate(across(where(is.numeric), ~ round(.x, digits)))
}

desc_stats <- function(x) {
  n_nonmiss <- sum(!is.na(x))
  tibble(
    Mean = if (n_nonmiss > 0) mean(x, na.rm = TRUE) else NA_real_,
    Median = if (n_nonmiss > 0) median(x, na.rm = TRUE) else NA_real_,
    `Std. Dev.` = if (n_nonmiss > 1) sd(x, na.rm = TRUE) else NA_real_,
    Min = if (n_nonmiss > 0) min(x, na.rm = TRUE) else NA_real_,
    Max = if (n_nonmiss > 0) max(x, na.rm = TRUE) else NA_real_,
    N = n_nonmiss
  )
}

winsor_vec <- function(x, p_low = 0.01, p_high = 0.99, min_non_na = 1) {
  n_nonmiss <- sum(!is.na(x))
  if (n_nonmiss < min_non_na) return(x)
  q_low  <- quantile(x, p_low,  na.rm = TRUE, type = 7)
  q_high <- quantile(x, p_high, na.rm = TRUE, type = 7)
  pmin(pmax(x, q_low), q_high)
}

winsor5 <- function(x) {
  n_nonmiss <- sum(!is.na(x))
  if (n_nonmiss < 20) return(x)
  p5  <- quantile(x, 0.05, na.rm = TRUE, type = 7)
  p95 <- quantile(x, 0.95, na.rm = TRUE, type = 7)
  pmin(pmax(x, p5), p95)
}

ensure_columns <- function(df, cols, fill = NA_real_) {
  for (cl in cols) {
    if (!cl %in% names(df)) df[[cl]] <- fill
  }
  df
}

nw_mean_tstat <- function(x, lag = 12, min_n = 30) {
  x <- as.numeric(x[!is.na(x)])
  n <- length(x)
  
  if (n == 0) {
    return(c(mean = NA_real_, se = NA_real_, tstat = NA_real_, n = 0))
  }
  
  if (n < min_n) {
    return(c(mean = mean(x), se = NA_real_, tstat = NA_real_, n = n))
  }
  
  fit <- lm(x ~ 1)
  se  <- sqrt(NeweyWest(fit, lag = lag, prewhite = FALSE)[1, 1])
  c(mean = mean(x), se = se, tstat = unname(coef(fit)[1]) / se, n = n)
}

# ------------------------------------------------------------------------------
# BUG FM CORRIGÉ :
# on ne garde que les variables présentes dans la formule
# ------------------------------------------------------------------------------
run_fm <- function(formula_str, data, lag = 12, min_obs_month = 20) {
  form <- as.formula(formula_str)
  vars_needed <- all.vars(form)
  
  monthly_coefs <- data %>%
    group_by(month_id) %>%
    group_modify(~ {
      d <- .x %>%
        dplyr::select(all_of(vars_needed)) %>%
        na.omit()
      
      if (nrow(d) < min_obs_month) return(tibble())
      
      fit <- tryCatch(lm(form, data = d), error = function(e) NULL)
      if (is.null(fit)) return(tibble())
      
      tibble(
        variable = names(coef(fit)),
        coef = as.numeric(coef(fit))
      )
    }) %>%
    ungroup() %>%
    pivot_wider(names_from = variable, values_from = coef)
  
  result_list <- lapply(setdiff(names(monthly_coefs), "month_id"), function(v) {
    s <- nw_mean_tstat(monthly_coefs[[v]], lag = lag)
    tibble(
      Variable = v,
      Coefficient = round(s["mean"], 6),
      tstat_NW = round(s["tstat"], 4),
      N_months = s["n"]
    )
  })
  
  bind_rows(result_list)
}

format_cell <- function(coef, tstat) {
  if (is.na(coef) || is.na(tstat)) return("—")
  stars <- dplyr::case_when(
    abs(tstat) >= 2.576 ~ "***",
    abs(tstat) >= 1.960 ~ "**",
    abs(tstat) >= 1.645 ~ "*",
    TRUE ~ ""
  )
  sprintf("%.4f (%.2f)%s", coef, tstat, stars)
}

make_academic_table <- function(fm_summary, variable_map, row_order) {
  fm_formatted <- fm_summary %>%
    mutate(
      Cell = mapply(format_cell, Coefficient, tstat_NW),
      Variable = dplyr::recode(Variable, !!!as.list(variable_map), .default = Variable)
    )
  
  n_months_row <- fm_formatted %>%
    group_by(Model) %>%
    summarise(N = first(N_months), .groups = "drop") %>%
    mutate(Variable = "N (months)", Cell = as.character(N)) %>%
    select(Model, Variable, Cell)
  
  bind_rows(
    fm_formatted %>% select(Model, Variable, Cell),
    n_months_row
  ) %>%
    pivot_wider(names_from = Model, values_from = Cell, values_fill = "—") %>%
    mutate(Variable = factor(Variable, levels = row_order)) %>%
    arrange(Variable)
}

build_ff_factors <- function(data) {
  factor_input <- data %>%
    group_by(month_id) %>%
    mutate(
      size_median = median(size_lag1, na.rm = TRUE),
      bm_p30 = quantile(bm_lag1, 0.30, na.rm = TRUE, type = 7),
      bm_p70 = quantile(bm_lag1, 0.70, na.rm = TRUE, type = 7),
      size_grp = if_else(size_lag1 <= size_median, "S", "B"),
      bm_grp = case_when(
        bm_lag1 <= bm_p30 ~ "L",
        bm_lag1 >  bm_p70 ~ "H",
        TRUE ~ "M"
      ),
      sb_grp = paste(size_grp, bm_grp, sep = "/")
    ) %>%
    ungroup()
  
  ff <- factor_input %>%
    filter(!is.na(sb_grp)) %>%
    group_by(month_id, sb_grp) %>%
    summarise(port_ret = mean(ret_lead1, na.rm = TRUE), .groups = "drop") %>%
    pivot_wider(names_from = sb_grp, values_from = port_ret) %>%
    arrange(month_id)
  
  ff <- ensure_columns(ff, c("S/L", "S/M", "S/H", "B/L", "B/M", "B/H"))
  
  ff %>%
    mutate(
      SMB = ((`S/L` + `S/M` + `S/H`) / 3) - ((`B/L` + `B/M` + `B/H`) / 3),
      HML = ((`S/H` + `B/H`) / 2) - ((`S/L` + `B/L`) / 2)
    )
}

get_alpha <- function(y, X_df, lag = 12) {
  fit <- lm(y ~ ., data = X_df)
  vc  <- NeweyWest(fit, lag = lag, prewhite = FALSE)
  coefs  <- coef(fit)
  ses    <- sqrt(diag(vc))
  tstats <- coefs / ses
  
  list(
    alpha = coefs[1],
    alpha_tstat = tstats[1],
    coefs = coefs,
    ses = ses,
    tstats = tstats
  )
}

build_market_states <- function(data) {
  data %>%
    group_by(month_id) %>%
    summarise(
      mkt_ret = mean(mkt_ret, na.rm = TRUE),
      rf = mean(rf, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    arrange(month_id) %>%
    mutate(
      mkt_excess = mkt_ret - rf,
      mkt_ret_lag = lag(mkt_ret),
      mkt_12m_past = slide_dbl(
        mkt_ret_lag,
        ~ if (any(is.na(.x))) NA_real_ else prod(1 + .x) - 1,
        .before = 11,
        .complete = TRUE
      ),
      market_state = case_when(
        is.na(mkt_12m_past) ~ NA_character_,
        mkt_12m_past > 0 ~ "Bull",
        mkt_12m_past <= 0 ~ "Bear",
        TRUE ~ NA_character_
      ),
      bear_dummy = case_when(
        is.na(mkt_12m_past) ~ NA_integer_,
        mkt_12m_past <= 0 ~ 1L,
        TRUE ~ 0L
      )
    ) %>%
    select(month_id, rf, mkt_ret, mkt_excess, mkt_12m_past, market_state, bear_dummy)
}

write_xlsx_safely <- function(sheets, path) {
  write_xlsx(lapply(sheets, as.data.frame), path = path)
}

# ==============================================================================
# 3. LOAD DATA
# ==============================================================================
monthly <- read_dta(PATH_MONTHLY)
daily   <- read_dta(PATH_DAILY)

cat("\n[INFO] monthly dimensions:\n"); print(dim(monthly))
cat("\n[INFO] daily dimensions:\n"); print(dim(daily))

# ==============================================================================
# 4. BUILD BASELINE PANEL
# ==============================================================================

# ------------------------------------------------------------------------------
# 4.1 Base monthly panel
# ------------------------------------------------------------------------------
panel_monthly_base <- monthly %>%
  transmute(
    date,
    month_id,
    dscd = as.character(dscd),
    stock_id,
    ret = ret_monthly_local_w,
    rf = rf_local,
    ret_excess = ret_monthly_local_w - rf_local,
    size = mv_monthly_local,
    log_size = log(mv_monthly_local),
    mtb = mtb_monthly,
    bm = 1 / mtb_monthly,
    mkt_ret = mkt_local_vw,
    index_ret = index_vw_local,
    price = prc_monthly_local,
    penny_stock,
    exname,
    ecname
  ) %>%
  filter(
    !is.na(dscd),
    !is.na(date),
    !is.na(ret),
    !is.na(size), size > 0,
    !is.na(mtb), mtb > 0,
    penny_stock != 1
  ) %>%
  arrange(dscd, month_id) %>%
  group_by(dscd) %>%
  mutate(
    ret_lead1 = lead(ret),
    ret_excess_lead1 = lead(ret_excess),
    size_lag1 = lag(size),
    log_size_lag1 = lag(log_size),
    bm_lag1 = lag(bm)
  ) %>%
  ungroup()

cat("\n[INFO] panel_monthly_base:\n"); print(dim(panel_monthly_base))

# ------------------------------------------------------------------------------
# 4.2 Monthly liquidity proxies from daily data
# ------------------------------------------------------------------------------
liq_monthly_from_daily <- daily %>%
  filter(penny_stock != 1) %>%
  transmute(
    dscd = as.character(dscd),
    month_id,
    ret_d = ret_daily_local_w,
    ret_d_raw = ret_daily_local,
    dollar_vol_d = dollar_volume_daily_w,
    turnover_d = turnover_daily_w,
    high_d = prc_high_daily_local,
    low_d = prc_low_daily_local
  ) %>%
  mutate(
    amihud_d = case_when(
      !is.na(ret_d) & !is.na(dollar_vol_d) & dollar_vol_d > 0 ~ abs(ret_d) / dollar_vol_d,
      TRUE ~ NA_real_
    ),
    zero_return_d = case_when(
      is.na(ret_d_raw) ~ NA_real_,
      ret_d_raw == 0 ~ 1,
      TRUE ~ 0
    ),
    hl_range_d = case_when(
      !is.na(high_d) & !is.na(low_d) & (high_d + low_d) > 0 ~
        (high_d - low_d) / ((high_d + low_d) / 2),
      TRUE ~ NA_real_
    )
  ) %>%
  group_by(dscd, month_id) %>%
  summarise(
    n_amihud_days   = sum(!is.na(amihud_d)),
    n_zero_days     = sum(!is.na(zero_return_d)),
    n_hl_days       = sum(!is.na(hl_range_d)),
    n_turnover_days = sum(!is.na(turnover_d)),
    amihud_m      = if_else(n_amihud_days   >= 15, mean(amihud_d, na.rm = TRUE), NA_real_),
    zero_return_m = if_else(n_zero_days     >= 15, mean(zero_return_d, na.rm = TRUE), NA_real_),
    hl_range_m    = if_else(n_hl_days       >= 15, mean(hl_range_d, na.rm = TRUE), NA_real_),
    turnover_m    = if_else(n_turnover_days >= 15, mean(turnover_d, na.rm = TRUE), NA_real_),
    .groups = "drop"
  )

cat("\n[INFO] liq_monthly_from_daily:\n"); print(dim(liq_monthly_from_daily))

# ------------------------------------------------------------------------------
# 4.3 Merge + baseline winsorisation (Amihud only)
# ------------------------------------------------------------------------------
panel_main <- panel_monthly_base %>%
  left_join(liq_monthly_from_daily, by = c("dscd", "month_id")) %>%
  arrange(dscd, month_id)

panel_main <- panel_main %>%
  group_by(month_id) %>%
  mutate(
    amihud_m_w1 = winsor_vec(amihud_m, p_low = 0.01, p_high = 0.99, min_non_na = 1)
  ) %>%
  ungroup() %>%
  arrange(dscd, month_id) %>%
  group_by(dscd) %>%
  mutate(
    amihud_lag1      = lag(amihud_m_w1),
    zero_return_lag1 = lag(zero_return_m),
    hl_range_lag1    = lag(hl_range_m),
    turnover_lag1    = lag(turnover_m)
  ) %>%
  ungroup()

cat("\n[INFO] panel_main:\n"); print(dim(panel_main))

# ------------------------------------------------------------------------------
# 4.4 Baseline sample
# ------------------------------------------------------------------------------
baseline_sample <- panel_main %>%
  filter(
    !is.na(ret_lead1),
    !is.na(ret_excess_lead1),
    !is.na(size_lag1), size_lag1 > 0,
    !is.na(log_size_lag1),
    !is.na(bm_lag1), bm_lag1 > 0,
    !is.na(amihud_lag1)
  )

cat("\n[INFO] baseline_sample:\n"); print(dim(baseline_sample))

baseline_sample %>%
  summarise(
    n_obs = n(),
    n_firms = n_distinct(dscd),
    n_months = n_distinct(month_id),
    first_month = min(month_id, na.rm = TRUE),
    last_month = max(month_id, na.rm = TRUE)
  ) %>% print()

saveRDS(baseline_sample, PATH_BASELINE_RDS)

# ==============================================================================
# 5. TABLES 1 À 6 — DESCRIPTIFS BASELINE
# ==============================================================================
panel_liq_desc <- panel_main %>%
  group_by(month_id) %>%
  mutate(
    amihud_w = winsor_vec(amihud_m, p_low = 0.01, p_high = 0.99, min_non_na = 1)
  ) %>%
  ungroup()

table_raw_liq <- bind_rows(
  desc_stats(liq_monthly_from_daily$amihud_m)      %>% mutate(Proxy = "Amihud"),
  desc_stats(liq_monthly_from_daily$zero_return_m) %>% mutate(Proxy = "Zero-return"),
  desc_stats(liq_monthly_from_daily$hl_range_m)    %>% mutate(Proxy = "High-low range"),
  desc_stats(liq_monthly_from_daily$turnover_m)    %>% mutate(Proxy = "Turnover")
) %>%
  select(Proxy, Mean, Median, `Std. Dev.`, Min, Max, N) %>%
  round_numeric()

table_before_winsor <- bind_rows(
  desc_stats(panel_main$amihud_m)      %>% mutate(Proxy = "Amihud"),
  desc_stats(panel_main$zero_return_m) %>% mutate(Proxy = "Zero-return"),
  desc_stats(panel_main$hl_range_m)    %>% mutate(Proxy = "High-low range"),
  desc_stats(panel_main$turnover_m)    %>% mutate(Proxy = "Turnover")
) %>%
  select(Proxy, Mean, Median, `Std. Dev.`, Min, Max, N) %>%
  round_numeric()

table_after_winsor <- bind_rows(
  desc_stats(panel_liq_desc$amihud_w)      %>% mutate(Proxy = "Amihud"),
  desc_stats(panel_liq_desc$zero_return_m) %>% mutate(Proxy = "Zero-return"),
  desc_stats(panel_liq_desc$hl_range_m)    %>% mutate(Proxy = "High-low range"),
  desc_stats(panel_liq_desc$turnover_m)    %>% mutate(Proxy = "Turnover")
) %>%
  select(Proxy, Mean, Median, `Std. Dev.`, Min, Max, N) %>%
  round_numeric()

panel_liq_desc <- panel_liq_desc %>%
  arrange(dscd, month_id) %>%
  group_by(dscd) %>%
  mutate(
    amihud_3m      = slide_dbl(amihud_w,      ~ mean(.x, na.rm = TRUE), .before = 2, .complete = TRUE),
    zero_return_3m = slide_dbl(zero_return_m, ~ mean(.x, na.rm = TRUE), .before = 2, .complete = TRUE),
    hl_range_3m    = slide_dbl(hl_range_m,    ~ mean(.x, na.rm = TRUE), .before = 2, .complete = TRUE),
    turnover_3m    = slide_dbl(turnover_m,    ~ mean(.x, na.rm = TRUE), .before = 2, .complete = TRUE)
  ) %>%
  ungroup()

table_3m <- bind_rows(
  desc_stats(panel_liq_desc$amihud_3m)      %>% mutate(Proxy = "Amihud (3-month MA)"),
  desc_stats(panel_liq_desc$zero_return_3m) %>% mutate(Proxy = "Zero-return (3-month MA)"),
  desc_stats(panel_liq_desc$hl_range_3m)    %>% mutate(Proxy = "High-low range (3-month MA)"),
  desc_stats(panel_liq_desc$turnover_3m)    %>% mutate(Proxy = "Turnover (3-month MA)")
) %>%
  select(Proxy, Mean, Median, `Std. Dev.`, Min, Max, N) %>%
  round_numeric()

table_descriptive_baseline <- bind_rows(
  ret_lead1        = desc_stats(baseline_sample$ret_lead1),
  ret_excess_lead1 = desc_stats(baseline_sample$ret_excess_lead1),
  size_lag1        = desc_stats(baseline_sample$size_lag1),
  log_size_lag1    = desc_stats(baseline_sample$log_size_lag1),
  bm_lag1          = desc_stats(baseline_sample$bm_lag1),
  amihud_lag1      = desc_stats(baseline_sample$amihud_lag1),
  zero_return_lag1 = desc_stats(baseline_sample$zero_return_lag1),
  hl_range_lag1    = desc_stats(baseline_sample$hl_range_lag1),
  turnover_lag1    = desc_stats(baseline_sample$turnover_lag1),
  .id = "Variable"
) %>%
  round_numeric()

size_comparison_sample <- baseline_sample %>%
  group_by(month_id) %>%
  mutate(
    size_median_month = median(size_lag1, na.rm = TRUE),
    size_group = case_when(
      size_lag1 <= size_median_month ~ "Small",
      size_lag1 > size_median_month ~ "Big",
      TRUE ~ NA_character_
    )
  ) %>%
  ungroup()

table_liquidity_by_size <- size_comparison_sample %>%
  group_by(size_group) %>%
  summarise(
    mean_amihud        = mean(amihud_lag1, na.rm = TRUE),
    median_amihud      = median(amihud_lag1, na.rm = TRUE),
    mean_zero_return   = mean(zero_return_lag1, na.rm = TRUE),
    median_zero_return = median(zero_return_lag1, na.rm = TRUE),
    mean_hl_range      = mean(hl_range_lag1, na.rm = TRUE),
    median_hl_range    = median(hl_range_lag1, na.rm = TRUE),
    mean_turnover      = mean(turnover_lag1, na.rm = TRUE),
    median_turnover    = median(turnover_lag1, na.rm = TRUE),
    n_obs = n(),
    .groups = "drop"
  ) %>%
  round_numeric()

write_xlsx_safely(
  list(
    "Table1_Raw_Construction" = table_raw_liq,
    "Table2_Before_Winsor" = table_before_winsor,
    "Table3_After_Winsor" = table_after_winsor,
    "Table4_Moving_Averages_3M" = table_3m,
    "Table5_Baseline_Sample" = table_descriptive_baseline,
    "Table6_Liquidity_by_Size" = table_liquidity_by_size
  ),
  file.path(OUT_DIR, "Tables_1_to_6_Baseline.xlsx")
)

# ==============================================================================
# 6. TABLES 7 À 9 — TRI UNIVARIÉ AMIHUD + CAPM / FF3
# ==============================================================================
portfolios_amihud <- baseline_sample %>%
  group_by(month_id) %>%
  mutate(
    q20 = quantile(amihud_lag1, 0.20, na.rm = TRUE, type = 7),
    q40 = quantile(amihud_lag1, 0.40, na.rm = TRUE, type = 7),
    q60 = quantile(amihud_lag1, 0.60, na.rm = TRUE, type = 7),
    q80 = quantile(amihud_lag1, 0.80, na.rm = TRUE, type = 7),
    quintile = case_when(
      amihud_lag1 <= q20 ~ 1,
      amihud_lag1 <= q40 ~ 2,
      amihud_lag1 <= q60 ~ 3,
      amihud_lag1 <= q80 ~ 4,
      amihud_lag1 >  q80 ~ 5,
      TRUE ~ NA_real_
    )
  ) %>%
  ungroup()

portfolio_returns <- portfolios_amihud %>%
  filter(!is.na(quintile)) %>%
  group_by(month_id, quintile) %>%
  summarise(
    port_ret = mean(ret_lead1, na.rm = TRUE),
    port_ret_excess = mean(ret_excess_lead1, na.rm = TRUE),
    rf_avg = mean(rf, na.rm = TRUE),
    mkt_ret_avg = mean(mkt_ret, na.rm = TRUE),
    n_stocks = n(),
    .groups = "drop"
  )

portfolio_returns_wide_excess <- portfolio_returns %>%
  select(month_id, quintile, port_ret_excess) %>%
  pivot_wider(names_from = quintile, values_from = port_ret_excess, names_prefix = "Q") %>%
  arrange(month_id)

portfolio_returns_wide_excess <- ensure_columns(portfolio_returns_wide_excess, paste0("Q", 1:5)) %>%
  mutate(Q5_minus_Q1 = Q5 - Q1)

factors_monthly <- baseline_sample %>%
  group_by(month_id) %>%
  summarise(
    rf = mean(rf, na.rm = TRUE),
    mkt_ret = mean(mkt_ret, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(mkt_excess = mkt_ret - rf) %>%
  arrange(month_id)

ff_portfolios <- build_ff_factors(baseline_sample)

ts_data <- portfolio_returns_wide_excess %>%
  left_join(factors_monthly, by = "month_id") %>%
  left_join(ff_portfolios %>% select(month_id, SMB, HML), by = "month_id") %>%
  filter(!is.na(mkt_excess), !is.na(SMB), !is.na(HML))

results_excess <- bind_rows(
  Q1 = nw_mean_tstat(ts_data$Q1),
  Q2 = nw_mean_tstat(ts_data$Q2),
  Q3 = nw_mean_tstat(ts_data$Q3),
  Q4 = nw_mean_tstat(ts_data$Q4),
  Q5 = nw_mean_tstat(ts_data$Q5),
  `Q5 - Q1` = nw_mean_tstat(ts_data$Q5_minus_Q1),
  .id = "Portfolio"
) %>%
  round_numeric()

capm_results <- list()
for (port in c("Q1", "Q2", "Q3", "Q4", "Q5", "Q5_minus_Q1")) {
  out <- get_alpha(ts_data[[port]], ts_data %>% select(mkt_excess))
  capm_results[[port]] <- tibble(
    Portfolio = port,
    alpha = out$alpha,
    alpha_tstat = out$alpha_tstat,
    beta_mkt = out$coefs[2],
    beta_mkt_tstat = out$tstats[2]
  )
}
capm_table <- bind_rows(capm_results) %>% round_numeric()

ff3_results <- list()
for (port in c("Q1", "Q2", "Q3", "Q4", "Q5", "Q5_minus_Q1")) {
  out <- get_alpha(ts_data[[port]], ts_data %>% select(mkt_excess, SMB, HML))
  ff3_results[[port]] <- tibble(
    Portfolio = port,
    alpha = out$alpha,
    alpha_tstat = out$alpha_tstat,
    beta_mkt = out$coefs[2],
    beta_smb = out$coefs[3],
    beta_hml = out$coefs[4],
    beta_mkt_tstat = out$tstats[2],
    beta_smb_tstat = out$tstats[3],
    beta_hml_tstat = out$tstats[4]
  )
}
ff3_table <- bind_rows(ff3_results) %>% round_numeric()

write_xlsx_safely(
  list(
    "Table7_MeanExcessReturns" = results_excess,
    "Table8_CAPM_Alphas" = capm_table,
    "Table9_FF3_Alphas" = ff3_table,
    "Portfolio_TS" = ts_data,
    "FF_Factors_TS" = ff_portfolios
  ),
  file.path(OUT_DIR, "Tables_7_to_9_Univariate_Amihud.xlsx")
)

# ==============================================================================
# 7. TABLES 10 À 13 — BIVARIATE SEQUENTIAL SORTS Size × Amihud (5x5)
# ==============================================================================
biv_sample <- baseline_sample %>%
  group_by(month_id) %>%
  mutate(
    s_q20 = quantile(size_lag1, 0.20, na.rm = TRUE, type = 7),
    s_q40 = quantile(size_lag1, 0.40, na.rm = TRUE, type = 7),
    s_q60 = quantile(size_lag1, 0.60, na.rm = TRUE, type = 7),
    s_q80 = quantile(size_lag1, 0.80, na.rm = TRUE, type = 7),
    size_quintile = case_when(
      size_lag1 <= s_q20 ~ 1,
      size_lag1 <= s_q40 ~ 2,
      size_lag1 <= s_q60 ~ 3,
      size_lag1 <= s_q80 ~ 4,
      size_lag1 >  s_q80 ~ 5,
      TRUE ~ NA_real_
    )
  ) %>%
  ungroup() %>%
  group_by(month_id, size_quintile) %>%
  mutate(
    a_q20 = quantile(amihud_lag1, 0.20, na.rm = TRUE, type = 7),
    a_q40 = quantile(amihud_lag1, 0.40, na.rm = TRUE, type = 7),
    a_q60 = quantile(amihud_lag1, 0.60, na.rm = TRUE, type = 7),
    a_q80 = quantile(amihud_lag1, 0.80, na.rm = TRUE, type = 7),
    amihud_quintile = case_when(
      amihud_lag1 <= a_q20 ~ 1,
      amihud_lag1 <= a_q40 ~ 2,
      amihud_lag1 <= a_q60 ~ 3,
      amihud_lag1 <= a_q80 ~ 4,
      amihud_lag1 >  a_q80 ~ 5,
      TRUE ~ NA_real_
    )
  ) %>%
  ungroup()

port_25 <- biv_sample %>%
  filter(!is.na(size_quintile), !is.na(amihud_quintile)) %>%
  group_by(month_id, size_quintile, amihud_quintile) %>%
  summarise(
    port_ret = mean(ret_lead1, na.rm = TRUE),
    port_ret_excess = mean(ret_excess_lead1, na.rm = TRUE),
    n_stocks = n(),
    .groups = "drop"
  )

fill_check <- port_25 %>%
  group_by(size_quintile, amihud_quintile) %>%
  summarise(
    n_months = n(),
    avg_n_stocks = mean(n_stocks),
    .groups = "drop"
  ) %>%
  round_numeric(2)

mean_matrix <- port_25 %>%
  group_by(size_quintile, amihud_quintile) %>%
  summarise(mean_excess = mean(port_ret_excess, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = amihud_quintile, values_from = mean_excess, names_prefix = "A") %>%
  round_numeric()

spread_within_size <- port_25 %>%
  select(month_id, size_quintile, amihud_quintile, port_ret_excess) %>%
  pivot_wider(names_from = amihud_quintile, values_from = port_ret_excess, names_prefix = "A") %>%
  arrange(month_id)

spread_within_size <- ensure_columns(spread_within_size, paste0("A", 1:5)) %>%
  mutate(A5_minus_A1 = A5 - A1)

spread_results_fixed <- spread_within_size %>%
  group_by(size_quintile) %>%
  group_modify(~ {
    s <- nw_mean_tstat(.x$A5_minus_A1, lag = 12)
    tibble(
      mean_A1 = round(mean(.x$A1, na.rm = TRUE), 6),
      mean_A5 = round(mean(.x$A5, na.rm = TRUE), 6),
      mean_spread = round(s["mean"], 6),
      tstat_NW = round(s["tstat"], 4),
      n_months = s["n"]
    )
  }) %>%
  ungroup()

size_neutral_spread <- spread_within_size %>%
  group_by(month_id) %>%
  summarise(avg_spread = mean(A5_minus_A1, na.rm = TRUE), .groups = "drop")

sn_summary <- nw_mean_tstat(size_neutral_spread$avg_spread, lag = 12)

sn_data <- size_neutral_spread %>%
  left_join(factors_monthly, by = "month_id") %>%
  left_join(ff_portfolios %>% select(month_id, SMB, HML), by = "month_id") %>%
  filter(!is.na(mkt_excess), !is.na(SMB), !is.na(HML))

capm_sn <- get_alpha(sn_data$avg_spread, sn_data %>% select(mkt_excess))
ff3_sn  <- get_alpha(sn_data$avg_spread, sn_data %>% select(mkt_excess, SMB, HML))

sn_alpha_table <- tibble(
  Model = c("CAPM", "Fama-French 3F"),
  alpha = c(capm_sn$coefs[1], ff3_sn$coefs[1]),
  alpha_tstat = c(capm_sn$tstats[1], ff3_sn$tstats[1]),
  beta_mkt = c(capm_sn$coefs[2], ff3_sn$coefs[2]),
  beta_mkt_tstat = c(capm_sn$tstats[2], ff3_sn$tstats[2]),
  beta_smb = c(NA, ff3_sn$coefs[3]),
  beta_smb_tstat = c(NA, ff3_sn$tstats[3]),
  beta_hml = c(NA, ff3_sn$coefs[4]),
  beta_hml_tstat = c(NA, ff3_sn$tstats[4])
) %>%
  round_numeric()

table_12_size_neutral <- tibble(
  Mean = round(sn_summary["mean"], 6),
  tstat_NW = round(sn_summary["tstat"], 4),
  N_months = sn_summary["n"]
)

write_xlsx_safely(
  list(
    "Table10_Mean_Returns_5x5" = mean_matrix,
    "Table11_Spread_by_Size" = spread_results_fixed,
    "Table12_SizeNeutral_Premium" = table_12_size_neutral,
    "Table12_SizeNeutral_TS" = size_neutral_spread,
    "Table13_SN_Alphas" = sn_alpha_table,
    "Fill_Check" = fill_check
  ),
  file.path(OUT_DIR, "Tables_10_to_13_Bivariate_Amihud.xlsx")
)

# ==============================================================================
# 8. TABLES 14 À 16 — BULL VS BEAR MARKETS
# ==============================================================================
market_ts <- build_market_states(baseline_sample)

portfolio_wide <- portfolio_returns %>%
  select(month_id, quintile, port_ret_excess) %>%
  pivot_wider(names_from = quintile, values_from = port_ret_excess, names_prefix = "Q") %>%
  arrange(month_id)

portfolio_wide <- ensure_columns(portfolio_wide, paste0("Q", 1:5)) %>%
  mutate(Q5_minus_Q1 = Q5 - Q1)

ts_state <- portfolio_wide %>%
  left_join(market_ts %>% select(month_id, market_state, mkt_excess), by = "month_id") %>%
  filter(!is.na(market_state))

table_state <- ts_state %>%
  select(market_state, Q1, Q2, Q3, Q4, Q5, Q5_minus_Q1) %>%
  pivot_longer(cols = -market_state, names_to = "Portfolio", values_to = "ret") %>%
  group_by(market_state, Portfolio) %>%
  group_modify(~ {
    s <- nw_mean_tstat(.x$ret, lag = 12)
    tibble(
      mean_ret = round(s["mean"], 6),
      tstat_NW = round(s["tstat"], 4),
      n_months = s["n"]
    )
  }) %>%
  ungroup()

stat_bull <- nw_mean_tstat(ts_state %>% filter(market_state == "Bull") %>% pull(Q5_minus_Q1), lag = 12)
stat_bear <- nw_mean_tstat(ts_state %>% filter(market_state == "Bear") %>% pull(Q5_minus_Q1), lag = 12)

ts_state_diff <- ts_state %>%
  mutate(bear_dummy = if_else(market_state == "Bear", 1, 0))

diff_fit <- lm(Q5_minus_Q1 ~ bear_dummy, data = ts_state_diff)
diff_vc  <- NeweyWest(diff_fit, lag = 12, prewhite = FALSE)
diff_tstats <- coef(diff_fit) / sqrt(diag(diff_vc))

table_diff <- tibble(
  State = c("Bull", "Bear", "Bear - Bull (diff)"),
  Mean_Spread = c(
    round(stat_bull["mean"], 6),
    round(stat_bear["mean"], 6),
    round(coef(diff_fit)[2], 6)
  ),
  tstat_NW = c(
    round(stat_bull["tstat"], 4),
    round(stat_bear["tstat"], 4),
    round(diff_tstats[2], 4)
  ),
  N_months = c(
    stat_bull["n"],
    stat_bear["n"],
    stat_bull["n"] + stat_bear["n"]
  )
)

ts_full <- ts_state %>%
  left_join(ff_portfolios %>% select(month_id, SMB, HML), by = "month_id") %>%
  filter(!is.na(SMB), !is.na(HML))

run_alpha_state <- function(state_label) {
  d <- ts_full %>% filter(market_state == state_label)
  
  fit_capm <- lm(Q5_minus_Q1 ~ mkt_excess, data = d)
  vc_capm  <- NeweyWest(fit_capm, lag = 12, prewhite = FALSE)
  
  fit_ff3 <- lm(Q5_minus_Q1 ~ mkt_excess + SMB + HML, data = d)
  vc_ff3  <- NeweyWest(fit_ff3, lag = 12, prewhite = FALSE)
  
  tibble(
    State = state_label,
    CAPM_alpha = round(coef(fit_capm)[1], 6),
    CAPM_tstat = round(coef(fit_capm)[1] / sqrt(vc_capm[1, 1]), 4),
    FF3_alpha = round(coef(fit_ff3)[1], 6),
    FF3_tstat = round(coef(fit_ff3)[1] / sqrt(vc_ff3[1, 1]), 4),
    FF3_beta_smb = round(coef(fit_ff3)[3], 4),
    FF3_smb_tstat = round(coef(fit_ff3)[3] / sqrt(vc_ff3[3, 3]), 4),
    N = nrow(d)
  )
}

table_alphas_state <- bind_rows(
  run_alpha_state("Bull"),
  run_alpha_state("Bear")
)

write_xlsx_safely(
  list(
    "Table14_Returns_by_State" = table_state,
    "Table15_Bull_vs_Bear_Diff" = table_diff,
    "Table16_Alphas_by_State" = table_alphas_state,
    "MarketStates_TS" = market_ts,
    "Portfolio_State_TS" = ts_full
  ),
  file.path(OUT_DIR, "Tables_14_to_16_BullBear.xlsx")
)

# ==============================================================================
# 9. TABLE 17 — FAMA-MACBETH BASELINE
# ==============================================================================
fm_data <- baseline_sample %>%
  left_join(market_ts %>% select(month_id, bear_dummy), by = "month_id") %>%
  mutate(
    log_amihud_lag1 = log(amihud_lag1 + 1e-8)
  ) %>%
  select(
    month_id, dscd, ret_lead1,
    amihud_lag1, log_amihud_lag1,
    log_size_lag1, bm_lag1, bear_dummy
  ) %>%
  filter(
    !is.na(amihud_lag1),
    !is.na(log_size_lag1),
    !is.na(bm_lag1)
  )

fm1 <- run_fm("ret_lead1 ~ amihud_lag1", fm_data)
fm2 <- run_fm("ret_lead1 ~ amihud_lag1 + log_size_lag1", fm_data)
fm3 <- run_fm("ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1", fm_data)
fm4 <- run_fm("ret_lead1 ~ log_amihud_lag1 + log_size_lag1 + bm_lag1", fm_data)

fm_data_state <- fm_data %>% filter(!is.na(bear_dummy))

fm5_bull <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_state %>% filter(bear_dummy == 0)
)

fm5_bear <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_state %>% filter(bear_dummy == 1)
)

fm_summary <- bind_rows(
  fm1 %>% mutate(Model = "(1) Amihud only"),
  fm2 %>% mutate(Model = "(2) Amihud + Size"),
  fm3 %>% mutate(Model = "(3) Baseline FM"),
  fm4 %>% mutate(Model = "(4) log(Amihud)"),
  fm5_bull %>% mutate(Model = "(5a) Bull markets only"),
  fm5_bear %>% mutate(Model = "(5b) Bear markets only")
) %>%
  select(Model, Variable, Coefficient, tstat_NW, N_months)

table17_academic <- make_academic_table(
  fm_summary = fm_summary,
  variable_map = c(
    "(Intercept)" = "Intercept",
    "amihud_lag1" = "Amihud",
    "log_amihud_lag1" = "log(Amihud)",
    "log_size_lag1" = "log(Size)",
    "bm_lag1" = "B/M"
  ),
  row_order = c("Intercept", "Amihud", "log(Amihud)", "log(Size)", "B/M", "N (months)")
)

write_xlsx_safely(
  list(
    "Table17_FM_Summary" = fm_summary,
    "Table17_Academic" = table17_academic,
    "Model1_Amihud_Only" = fm1,
    "Model2_Amihud_Size" = fm2,
    "Model3_Baseline" = fm3,
    "Model4_LogAmihud" = fm4,
    "Model5a_Bull" = fm5_bull,
    "Model5b_Bear" = fm5_bear
  ),
  file.path(OUT_DIR, "Table_17_FamaMacBeth.xlsx")
)

cat("\n============================\n")
cat("BLOC 1/3 TERMINÉ\n")
cat("============================\n")
# ==============================================================================
# BLOC 2/3 — SAFE MANAGER SENTIMENT + TABLES 18 TO 21
# ==============================================================================

# Vérification baseline
if (!file.exists(PATH_BASELINE_RDS)) {
  stop("Le fichier baseline_sample_final.rds est introuvable. Exécute d'abord le BLOC 1.")
}

baseline_sample <- readRDS(PATH_BASELINE_RDS)

# ==============================================================================
# 10. SAFE MANAGER SENTIMENT
# ==============================================================================
sentiment_safe <- read_csv(PATH_SAFE, show_col_types = FALSE) %>%
  rename(sentiment = SAFE_Manager_Sentiment) %>%
  mutate(
    year_safe  = as.integer(sub("m.*", "", YYYYMM)),
    month_safe = as.integer(sub(".*m", "", YYYYMM)),
    date_safe  = as.Date(sprintf("%04d-%02d-01", year_safe, month_safe)),
    year_month = sprintf("%04d-%02d", year_safe, month_safe)
  ) %>%
  filter(!is.na(sentiment))

month_mapping <- monthly %>%
  transmute(
    month_id,
    year_month = format(as.Date(date), "%Y-%m")
  ) %>%
  distinct(month_id, year_month) %>%
  arrange(month_id)

sentiment_panel <- sentiment_safe %>%
  inner_join(month_mapping, by = "year_month") %>%
  arrange(month_id) %>%
  select(month_id, year_month, date_safe, sentiment)

if (nrow(sentiment_panel) == 0) {
  stop("Aucun match entre SAFE et month_id. Vérifie le parsing de year_month.")
}

sent_median_val <- median(sentiment_panel$sentiment, na.rm = TRUE)
sent_p33_val    <- quantile(sentiment_panel$sentiment, 0.33, na.rm = TRUE)
sent_p67_val    <- quantile(sentiment_panel$sentiment, 0.67, na.rm = TRUE)

sentiment_panel <- sentiment_panel %>%
  mutate(
    sent_state_2 = if_else(sentiment > sent_median_val, "High", "Low"),
    sent_state_3 = case_when(
      sentiment <= sent_p33_val ~ "Low",
      sentiment >= sent_p67_val ~ "High",
      TRUE ~ "Medium"
    ),
    high_sent_dummy = if_else(sent_state_2 == "High", 1L, 0L)
  )

baseline_sample_sent_only <- baseline_sample %>%
  left_join(
    sentiment_panel %>% select(month_id, sentiment, sent_state_2, sent_state_3, high_sent_dummy),
    by = "month_id"
  ) %>%
  filter(!is.na(sentiment), !is.na(sent_state_2))

saveRDS(baseline_sample_sent_only, PATH_BASELINE_SENT_RDS)

cat("\n[INFO] baseline_sample_sent_only:\n")
baseline_sample_sent_only %>%
  summarise(
    n_obs = n(),
    n_firms = n_distinct(dscd),
    n_months = n_distinct(month_id),
    first_month = min(month_id),
    last_month = max(month_id)
  ) %>% print()

# ------------------------------------------------------------------------------
# Table 18 — Mean excess returns by sentiment state
# ------------------------------------------------------------------------------
portfolios_sent <- baseline_sample_sent_only %>%
  group_by(month_id) %>%
  mutate(
    q20 = quantile(amihud_lag1, 0.20, na.rm = TRUE, type = 7),
    q40 = quantile(amihud_lag1, 0.40, na.rm = TRUE, type = 7),
    q60 = quantile(amihud_lag1, 0.60, na.rm = TRUE, type = 7),
    q80 = quantile(amihud_lag1, 0.80, na.rm = TRUE, type = 7),
    quintile = case_when(
      amihud_lag1 <= q20 ~ 1,
      amihud_lag1 <= q40 ~ 2,
      amihud_lag1 <= q60 ~ 3,
      amihud_lag1 <= q80 ~ 4,
      amihud_lag1 > q80 ~ 5,
      TRUE ~ NA_real_
    )
  ) %>%
  ungroup()

portfolio_returns_sent <- portfolios_sent %>%
  filter(!is.na(quintile)) %>%
  group_by(month_id, quintile) %>%
  summarise(
    port_ret_excess = mean(ret_excess_lead1, na.rm = TRUE),
    n_stocks = n(),
    .groups = "drop"
  )

port_wide_sent <- portfolio_returns_sent %>%
  select(month_id, quintile, port_ret_excess) %>%
  pivot_wider(names_from = quintile, values_from = port_ret_excess, names_prefix = "Q") %>%
  arrange(month_id)

port_wide_sent <- ensure_columns(port_wide_sent, paste0("Q", 1:5)) %>%
  mutate(Q5_minus_Q1 = Q5 - Q1)

ts_sent <- port_wide_sent %>%
  left_join(sentiment_panel %>% select(month_id, sent_state_2, sent_state_3), by = "month_id") %>%
  filter(!is.na(sent_state_2))

table_18 <- ts_sent %>%
  select(sent_state_2, Q1, Q2, Q3, Q4, Q5, Q5_minus_Q1) %>%
  pivot_longer(cols = -sent_state_2, names_to = "Portfolio", values_to = "ret") %>%
  group_by(sent_state_2, Portfolio) %>%
  group_modify(~ {
    s <- nw_mean_tstat(.x$ret, lag = 12)
    tibble(
      mean_ret = round(s["mean"], 6),
      tstat_NW = round(s["tstat"], 4),
      n_months = s["n"]
    )
  }) %>%
  ungroup() %>%
  mutate(
    Portfolio = factor(Portfolio, levels = c("Q1", "Q2", "Q3", "Q4", "Q5", "Q5_minus_Q1"))
  ) %>%
  arrange(sent_state_2, Portfolio)

# ------------------------------------------------------------------------------
# Table 19 — Spread High vs Low sentiment
# ------------------------------------------------------------------------------
stat_high <- nw_mean_tstat(
  ts_sent %>% filter(sent_state_2 == "High") %>% pull(Q5_minus_Q1),
  lag = 12
)

stat_low <- nw_mean_tstat(
  ts_sent %>% filter(sent_state_2 == "Low") %>% pull(Q5_minus_Q1),
  lag = 12
)

ts_sent_diff <- ts_sent %>%
  mutate(high_dummy = if_else(sent_state_2 == "High", 1, 0))

diff_fit <- lm(Q5_minus_Q1 ~ high_dummy, data = ts_sent_diff)
diff_vc  <- NeweyWest(diff_fit, lag = 12, prewhite = FALSE)
diff_tstats <- coef(diff_fit) / sqrt(diag(diff_vc))

table_19 <- tibble(
  State = c("Low Sentiment", "High Sentiment", "High - Low (diff)"),
  Mean_Spread = c(
    round(stat_low["mean"], 6),
    round(stat_high["mean"], 6),
    round(coef(diff_fit)[2], 6)
  ),
  tstat_NW = c(
    round(stat_low["tstat"], 4),
    round(stat_high["tstat"], 4),
    round(diff_tstats[2], 4)
  ),
  N_months = c(
    stat_low["n"],
    stat_high["n"],
    stat_low["n"] + stat_high["n"]
  )
)

# ------------------------------------------------------------------------------
# Table 20 — Conditional alphas by sentiment state
# ------------------------------------------------------------------------------
factors_monthly_sent <- baseline_sample_sent_only %>%
  group_by(month_id) %>%
  summarise(
    rf = mean(rf, na.rm = TRUE),
    mkt_ret = mean(mkt_ret, na.rm = TRUE),
    .groups = "drop"
  ) %>%
  mutate(mkt_excess = mkt_ret - rf)

ff_portfolios_sent <- build_ff_factors(baseline_sample_sent_only)

ts_full_sent <- ts_sent %>%
  left_join(factors_monthly_sent, by = "month_id") %>%
  left_join(ff_portfolios_sent %>% select(month_id, SMB, HML), by = "month_id") %>%
  filter(!is.na(mkt_excess), !is.na(SMB), !is.na(HML))

run_alpha_sent <- function(state_label) {
  d <- ts_full_sent %>% filter(sent_state_2 == state_label)
  
  fit_capm <- lm(Q5_minus_Q1 ~ mkt_excess, data = d)
  vc_capm  <- NeweyWest(fit_capm, lag = 12, prewhite = FALSE)
  
  fit_ff3 <- lm(Q5_minus_Q1 ~ mkt_excess + SMB + HML, data = d)
  vc_ff3  <- NeweyWest(fit_ff3, lag = 12, prewhite = FALSE)
  
  tibble(
    Sentiment_State = state_label,
    CAPM_alpha = round(coef(fit_capm)[1], 6),
    CAPM_tstat = round(coef(fit_capm)[1] / sqrt(vc_capm[1, 1]), 4),
    FF3_alpha = round(coef(fit_ff3)[1], 6),
    FF3_tstat = round(coef(fit_ff3)[1] / sqrt(vc_ff3[1, 1]), 4),
    FF3_beta_smb = round(coef(fit_ff3)[3], 4),
    FF3_smb_tstat = round(coef(fit_ff3)[3] / sqrt(vc_ff3[3, 3]), 4),
    FF3_beta_hml = round(coef(fit_ff3)[4], 4),
    FF3_hml_tstat = round(coef(fit_ff3)[4] / sqrt(vc_ff3[4, 4]), 4),
    N_months = nrow(d)
  )
}

table_20 <- bind_rows(
  run_alpha_sent("Low"),
  run_alpha_sent("High")
)

# ------------------------------------------------------------------------------
# Table 21 — Fama-MacBeth conditioned on sentiment
# ------------------------------------------------------------------------------
fm_data_sent <- baseline_sample_sent_only %>%
  mutate(
    log_amihud_lag1 = log(amihud_lag1 + 1e-8)
  ) %>%
  select(
    month_id, dscd, ret_lead1,
    amihud_lag1, log_amihud_lag1,
    log_size_lag1, bm_lag1,
    sentiment, high_sent_dummy
  ) %>%
  filter(
    !is.na(amihud_lag1),
    !is.na(log_size_lag1),
    !is.na(bm_lag1),
    !is.na(sentiment)
  )

fm_s1 <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_sent
)

fm_s2 <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_sent %>% filter(high_sent_dummy == 1)
)

fm_s3 <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_sent %>% filter(high_sent_dummy == 0)
)

fm_s4 <- run_fm(
  "ret_lead1 ~ log_amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_sent
)

fm_s5 <- run_fm(
  "ret_lead1 ~ log_amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_sent %>% filter(high_sent_dummy == 1)
)

fm_s6 <- run_fm(
  "ret_lead1 ~ log_amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_sent %>% filter(high_sent_dummy == 0)
)

table_21 <- bind_rows(
  fm_s1 %>% mutate(Model = "(S1) Baseline period sent"),
  fm_s2 %>% mutate(Model = "(S2) High Sentiment"),
  fm_s3 %>% mutate(Model = "(S3) Low Sentiment"),
  fm_s4 %>% mutate(Model = "(S4) log(Amihud) baseline"),
  fm_s5 %>% mutate(Model = "(S5) log(Amihud) High Sent"),
  fm_s6 %>% mutate(Model = "(S6) log(Amihud) Low Sent")
) %>%
  select(Model, Variable, Coefficient, tstat_NW, N_months)

table21_academic <- make_academic_table(
  fm_summary = table_21,
  variable_map = c(
    "(Intercept)" = "Intercept",
    "amihud_lag1" = "Amihud",
    "log_amihud_lag1" = "log(Amihud)",
    "log_size_lag1" = "log(Size)",
    "bm_lag1" = "B/M"
  ),
  row_order = c("Intercept", "Amihud", "log(Amihud)", "log(Size)", "B/M", "N (months)")
)

write_xlsx_safely(
  list(
    "Table18_Returns_by_Sentiment" = table_18,
    "Table19_Spread_HighLow" = table_19,
    "Table20_Alphas_by_Sentiment" = table_20,
    "Table21_FM_Sentiment" = table_21,
    "Table21_Academic" = table21_academic,
    "Sentiment_TS" = sentiment_panel,
    "Portfolio_Sent_TS" = ts_full_sent
  ),
  file.path(OUT_DIR, "Tables_18_to_21_Sentiment.xlsx")
)

cat("\n============================\n")
cat("BLOC 2/3 TERMINÉ\n")
cat("============================\n")
# ==============================================================================
# BLOC 3/3 — EDGE + TABLE 22 + ROBUSTNESS 5%
# ==============================================================================

# Vérification baseline
if (!file.exists(PATH_BASELINE_RDS)) {
  stop("Le fichier baseline_sample_final.rds est introuvable. Exécute d'abord le BLOC 1.")
}

baseline_sample <- readRDS(PATH_BASELINE_RDS)

# ==============================================================================
# 11. EDGE PROXY — CONSTRUCTION + DESCRIPTIFS + CORRÉLATIONS + FM + SORTS
# ==============================================================================
if (!require("bidask")) install.packages("bidask")
library(bidask)

cat("\n=== DIAGNOSTIC : colonnes prix dans daily ===\n")
prc_cols <- grep("prc|open|close|high|low", names(daily), value = TRUE, ignore.case = TRUE)
print(prc_cols)

edge_input <- daily %>%
  filter(penny_stock != 1) %>%
  arrange(dscd, date) %>%
  group_by(dscd) %>%
  mutate(
    open_proxy = lag(prc_daily_local)
  ) %>%
  ungroup() %>%
  transmute(
    dscd = as.character(dscd),
    month_id,
    date,
    open  = open_proxy,
    high  = prc_high_daily_local,
    low   = prc_low_daily_local,
    close = prc_daily_local
  ) %>%
  filter(!is.na(open), !is.na(high), !is.na(low), !is.na(close)) %>%
  arrange(dscd, date)

cat("\n[INFO] edge_input dimensions:\n"); print(dim(edge_input))

edge_monthly <- edge_input %>%
  group_by(dscd, month_id) %>%
  summarise(
    n_edge_days = n(),
    edge_m = if (n_edge_days >= 15) {
      tryCatch(
        bidask::edge(open = open, high = high, low = low, close = close, sign = FALSE),
        error = function(e) NA_real_
      )
    } else NA_real_,
    .groups = "drop"
  )

cat("\n[INFO] edge_monthly dimensions:\n"); print(dim(edge_monthly))
cat("\n[INFO] Distribution EDGE:\n"); print(summary(edge_monthly$edge_m))

panel_main_edge <- panel_main %>%
  left_join(edge_monthly %>% select(dscd, month_id, edge_m, n_edge_days), by = c("dscd", "month_id")) %>%
  group_by(month_id) %>%
  mutate(
    edge_m_w1 = winsor_vec(edge_m, p_low = 0.01, p_high = 0.99, min_non_na = 1)
  ) %>%
  ungroup() %>%
  arrange(dscd, month_id) %>%
  group_by(dscd) %>%
  mutate(
    edge_lag1 = lag(edge_m_w1)
  ) %>%
  ungroup()

baseline_sample_with_edge <- panel_main_edge %>%
  filter(
    !is.na(ret_lead1),
    !is.na(ret_excess_lead1),
    !is.na(size_lag1), size_lag1 > 0,
    !is.na(log_size_lag1),
    !is.na(bm_lag1), bm_lag1 > 0,
    !is.na(amihud_lag1)
  )

baseline_sample_edge_only <- baseline_sample_with_edge %>%
  filter(!is.na(edge_lag1))

saveRDS(baseline_sample_with_edge, PATH_BASELINE_EDGE_RDS)
saveRDS(baseline_sample_edge_only, PATH_BASELINE_EDGE_ONLYRDS)

table_edge_desc <- bind_rows(
  desc_stats(edge_monthly$edge_m) %>% mutate(Stage = "1. Raw construction"),
  desc_stats(panel_main_edge$edge_m) %>% mutate(Stage = "2. After merge (before winsor)"),
  desc_stats(panel_main_edge$edge_m_w1) %>% mutate(Stage = "3. After winsor 1%/99%")
) %>%
  select(Stage, Mean, Median, `Std. Dev.`, Min, Max, N) %>%
  round_numeric()

table_all_proxies <- bind_rows(
  desc_stats(panel_main_edge$amihud_m_w1)   %>% mutate(Proxy = "Amihud (winsor)"),
  desc_stats(panel_main_edge$zero_return_m) %>% mutate(Proxy = "Zero-return"),
  desc_stats(panel_main_edge$hl_range_m)    %>% mutate(Proxy = "High-low range"),
  desc_stats(panel_main_edge$turnover_m)    %>% mutate(Proxy = "Turnover"),
  desc_stats(panel_main_edge$edge_m_w1)     %>% mutate(Proxy = "EDGE (winsor)")
) %>%
  select(Proxy, Mean, Median, `Std. Dev.`, Min, Max, N) %>%
  round_numeric()

cor_data <- baseline_sample_with_edge %>%
  select(amihud_lag1, zero_return_lag1, hl_range_lag1, turnover_lag1, edge_lag1) %>%
  rename(
    Amihud = amihud_lag1,
    `Zero-ret` = zero_return_lag1,
    `HL range` = hl_range_lag1,
    Turnover = turnover_lag1,
    EDGE = edge_lag1
  )

cor_pearson  <- cor(cor_data, use = "pairwise.complete.obs", method = "pearson")
cor_spearman <- cor(cor_data, use = "pairwise.complete.obs", method = "spearman")

fm_data_edge <- baseline_sample_edge_only %>%
  mutate(
    log_edge_lag1 = log(edge_lag1 + 1e-8)
  ) %>%
  select(
    month_id, dscd, ret_lead1,
    edge_lag1, log_edge_lag1,
    log_size_lag1, bm_lag1
  ) %>%
  filter(!is.na(edge_lag1), !is.na(log_size_lag1), !is.na(bm_lag1))

fm_e1 <- run_fm("ret_lead1 ~ edge_lag1", fm_data_edge)
fm_e2 <- run_fm("ret_lead1 ~ edge_lag1 + log_size_lag1", fm_data_edge)
fm_e3 <- run_fm("ret_lead1 ~ edge_lag1 + log_size_lag1 + bm_lag1", fm_data_edge)
fm_e4 <- run_fm("ret_lead1 ~ log_edge_lag1 + log_size_lag1 + bm_lag1", fm_data_edge)

fm_edge_summary <- bind_rows(
  fm_e1 %>% mutate(Model = "(E1) EDGE only"),
  fm_e2 %>% mutate(Model = "(E2) EDGE + Size"),
  fm_e3 %>% mutate(Model = "(E3) Baseline EDGE"),
  fm_e4 %>% mutate(Model = "(E4) log(EDGE)")
) %>%
  select(Model, Variable, Coefficient, tstat_NW, N_months)

table_edge_academic <- make_academic_table(
  fm_summary = fm_edge_summary,
  variable_map = c(
    "(Intercept)" = "Intercept",
    "edge_lag1" = "EDGE",
    "log_edge_lag1" = "log(EDGE)",
    "log_size_lag1" = "log(Size)",
    "bm_lag1" = "B/M"
  ),
  row_order = c("Intercept", "EDGE", "log(EDGE)", "log(Size)", "B/M", "N (months)")
)

portfolios_edge <- baseline_sample_edge_only %>%
  group_by(month_id) %>%
  mutate(
    e_q20 = quantile(edge_lag1, 0.20, na.rm = TRUE, type = 7),
    e_q40 = quantile(edge_lag1, 0.40, na.rm = TRUE, type = 7),
    e_q60 = quantile(edge_lag1, 0.60, na.rm = TRUE, type = 7),
    e_q80 = quantile(edge_lag1, 0.80, na.rm = TRUE, type = 7),
    quintile_edge = case_when(
      edge_lag1 <= e_q20 ~ 1,
      edge_lag1 <= e_q40 ~ 2,
      edge_lag1 <= e_q60 ~ 3,
      edge_lag1 <= e_q80 ~ 4,
      edge_lag1 >  e_q80 ~ 5,
      TRUE ~ NA_real_
    )
  ) %>%
  ungroup()

portfolio_returns_edge <- portfolios_edge %>%
  filter(!is.na(quintile_edge)) %>%
  group_by(month_id, quintile_edge) %>%
  summarise(
    port_ret_excess = mean(ret_excess_lead1, na.rm = TRUE),
    n_stocks = n(),
    .groups = "drop"
  )

port_wide_edge <- portfolio_returns_edge %>%
  select(month_id, quintile_edge, port_ret_excess) %>%
  pivot_wider(names_from = quintile_edge, values_from = port_ret_excess, names_prefix = "Q") %>%
  arrange(month_id)

port_wide_edge <- ensure_columns(port_wide_edge, paste0("Q", 1:5)) %>%
  mutate(Q5_minus_Q1 = Q5 - Q1)

results_edge <- bind_rows(
  Q1 = nw_mean_tstat(port_wide_edge$Q1),
  Q2 = nw_mean_tstat(port_wide_edge$Q2),
  Q3 = nw_mean_tstat(port_wide_edge$Q3),
  Q4 = nw_mean_tstat(port_wide_edge$Q4),
  Q5 = nw_mean_tstat(port_wide_edge$Q5),
  `Q5 - Q1` = nw_mean_tstat(port_wide_edge$Q5_minus_Q1),
  .id = "Portfolio"
) %>%
  round_numeric()

write_xlsx_safely(
  list(
    "EDGE_Descriptives" = table_edge_desc,
    "All_Proxies_Desc" = table_all_proxies,
    "Correlation_Pearson" = as.data.frame(round(cor_pearson, 4)) %>% tibble::rownames_to_column("Proxy"),
    "Correlation_Spearman" = as.data.frame(round(cor_spearman, 4)) %>% tibble::rownames_to_column("Proxy"),
    "FM_EDGE_Summary" = fm_edge_summary,
    "FM_EDGE_Academic" = table_edge_academic,
    "EDGE_Univariate_Sorts" = results_edge,
    "Portfolio_TS_EDGE" = port_wide_edge
  ),
  file.path(OUT_DIR, "EDGE_Complete_Analysis.xlsx")
)

# ==============================================================================
# 12. TABLE 22 — ROBUSTESSE : FM BASELINE AVEC CHAQUE PROXY
# ==============================================================================
fm_data_robust <- baseline_sample %>%
  filter(!is.na(amihud_lag1), !is.na(log_size_lag1), !is.na(bm_lag1))

fm_r1 <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  fm_data_robust
)

fm_r2 <- run_fm(
  "ret_lead1 ~ zero_return_lag1 + log_size_lag1 + bm_lag1",
  fm_data_robust %>% filter(!is.na(zero_return_lag1))
)

fm_r3 <- run_fm(
  "ret_lead1 ~ hl_range_lag1 + log_size_lag1 + bm_lag1",
  fm_data_robust %>% filter(!is.na(hl_range_lag1))
)

fm_r4 <-
  fm_r4 <- run_fm(
    "ret_lead1 ~ turnover_lag1 + log_size_lag1 + bm_lag1",
    fm_data_robust %>% filter(!is.na(turnover_lag1))
  )

table_22 <- bind_rows(
  fm_r1 %>% mutate(Model = "(R1) Amihud (baseline)"),
  fm_r2 %>% mutate(Model = "(R2) Zero-return"),
  fm_r3 %>% mutate(Model = "(R3) High-low range"),
  fm_r4 %>% mutate(Model = "(R4) Turnover")
) %>%
  select(Model, Variable, Coefficient, tstat_NW, N_months)

table22_academic <- make_academic_table(
  fm_summary = table_22,
  variable_map = c(
    "(Intercept)" = "Intercept",
    "amihud_lag1" = "Liquidity proxy",
    "zero_return_lag1" = "Liquidity proxy",
    "hl_range_lag1" = "Liquidity proxy",
    "turnover_lag1" = "Liquidity proxy",
    "log_size_lag1" = "log(Size)",
    "bm_lag1" = "B/M"
  ),
  row_order = c("Intercept", "Liquidity proxy", "log(Size)", "B/M", "N (months)")
)

write_xlsx_safely(
  list(
    "Table22_Robustness_Proxies" = table_22,
    "Table22_Academic" = table22_academic
  ),
  file.path(OUT_DIR, "Table22_Robustness_Proxies.xlsx")
)

# ==============================================================================
# 13. ROBUSTNESS CHECK — 5% WINSORIZATION FOR ALL CONTINUOUS VARIABLES
# BUG DES LAGS CORRIGÉ :
# d'abord les lags par firme (dscd), ensuite winsorisation par mois
# ==============================================================================

# ------------------------------------------------------------------------------
# 13.1 Rebuild robust sample with correct lag structure
# ------------------------------------------------------------------------------
robust_panel <- panel_monthly_base %>%
  left_join(liq_monthly_from_daily, by = c("dscd", "month_id")) %>%
  arrange(dscd, month_id) %>%
  group_by(dscd) %>%
  mutate(
    amihud_lag1_raw      = lag(amihud_m),
    turnover_lag1_raw    = lag(turnover_m),
    hl_range_lag1_raw    = lag(hl_range_m),
    zero_return_lag1     = lag(zero_return_m)
  ) %>%
  ungroup() %>%
  group_by(month_id) %>%
  mutate(
    amihud_lag1   = winsor5(amihud_lag1_raw),
    turnover_lag1 = winsor5(turnover_lag1_raw),
    hl_range_lag1 = winsor5(hl_range_lag1_raw),
    log_size_lag1 = winsor5(log_size_lag1),
    bm_lag1       = winsor5(bm_lag1)
  ) %>%
  ungroup()

robust_sample <- robust_panel %>%
  filter(
    !is.na(ret_lead1),
    !is.na(amihud_lag1),
    !is.na(log_size_lag1),
    !is.na(bm_lag1)
  )

cat("\n[INFO] robust_sample dimensions:\n")
print(dim(robust_sample))

# ------------------------------------------------------------------------------
# 13.2 EDGE with 5% winsorization — corrected lag structure
# ------------------------------------------------------------------------------
edge_monthly_robust <- daily %>%
  filter(penny_stock != 1) %>%
  arrange(dscd, date) %>%
  group_by(dscd) %>%
  mutate(open_proxy = lag(prc_daily_local)) %>%
  ungroup() %>%
  transmute(
    dscd = as.character(dscd),
    month_id,
    date,
    open  = open_proxy,
    high  = prc_high_daily_local,
    low   = prc_low_daily_local,
    close = prc_daily_local
  ) %>%
  filter(!is.na(open), !is.na(high), !is.na(low), !is.na(close)) %>%
  group_by(dscd, month_id) %>%
  summarise(
    n_edge_days = n(),
    edge_m = if (n_edge_days >= 15) {
      tryCatch(
        bidask::edge(open = open, high = high, low = low, close = close, sign = FALSE),
        error = function(e) NA_real_
      )
    } else NA_real_,
    .groups = "drop"
  )

robust_sample <- robust_sample %>%
  left_join(edge_monthly_robust %>% select(dscd, month_id, edge_m), by = c("dscd", "month_id")) %>%
  arrange(dscd, month_id) %>%
  group_by(dscd) %>%
  mutate(
    edge_lag1_raw = lag(edge_m)
  ) %>%
  ungroup() %>%
  group_by(month_id) %>%
  mutate(
    edge_lag1 = winsor5(edge_lag1_raw)
  ) %>%
  ungroup()

# ------------------------------------------------------------------------------
# 13.3 Fama-MacBeth robustness regressions (5% winsorization)
# ------------------------------------------------------------------------------
fm_amihud_5pct <- run_fm(
  "ret_lead1 ~ amihud_lag1 + log_size_lag1 + bm_lag1",
  robust_sample
)

fm_edge_5pct <- run_fm(
  "ret_lead1 ~ edge_lag1 + log_size_lag1 + bm_lag1",
  robust_sample %>% filter(!is.na(edge_lag1))
)

table_robust_5pct <- bind_rows(
  fm_amihud_5pct %>% mutate(Model = "Amihud 5pct"),
  fm_edge_5pct %>% mutate(Model = "EDGE 5pct")
) %>%
  select(Model, Variable, Coefficient, tstat_NW, N_months)

table_robust_5pct_academic <- make_academic_table(
  fm_summary = table_robust_5pct,
  variable_map = c(
    "(Intercept)" = "Intercept",
    "amihud_lag1" = "Amihud",
    "edge_lag1" = "EDGE",
    "log_size_lag1" = "log(Size)",
    "bm_lag1" = "B/M"
  ),
  row_order = c("Intercept", "Amihud", "EDGE", "log(Size)", "B/M", "N (months)")
)

write_xlsx_safely(
  list(
    "Robustness_5pct_Summary" = table_robust_5pct,
    "Robustness_5pct_Academic" = table_robust_5pct_academic,
    "FM_Amihud_5pct" = fm_amihud_5pct,
    "FM_EDGE_5pct" = fm_edge_5pct
  ),
  file.path(OUT_DIR, "Robustesse_5pct_Results.xlsx")
)

cat("\n============================\n")
cat("BLOC 3/3 TERMINÉ\n")
cat("============================\n")
cat("\nFichiers créés :\n")
cat("- EDGE_Complete_Analysis.xlsx\n")
cat("- Table22_Robustness_Proxies.xlsx\n")
cat("- Robustesse_5pct_Results.xlsx\n")
############################################################
# BLOC FINAL : REGROUPER TOUS LES TABLEAUX DANS UN SEUL EXCEL
############################################################

# Packages
req_merge <- c("openxlsx", "readxl")
new_merge <- req_merge[!(req_merge %in% installed.packages()[, "Package"])]
if (length(new_merge) > 0) install.packages(new_merge)

library(openxlsx)
library(readxl)

OUT_DIR <- "C:/Users/najwa/Downloads"

# Fichiers Excel déjà générés par ton script
source_files <- c(
  file.path(OUT_DIR, "Tables_1_to_6_Baseline.xlsx"),
  file.path(OUT_DIR, "Tables_7_to_9_Univariate_Amihud.xlsx"),
  file.path(OUT_DIR, "Tables_10_to_13_Bivariate_Amihud.xlsx"),
  file.path(OUT_DIR, "Tables_14_to_16_BullBear.xlsx"),
  file.path(OUT_DIR, "Table_17_FamaMacBeth.xlsx"),
  file.path(OUT_DIR, "Tables_18_to_21_Sentiment.xlsx"),
  file.path(OUT_DIR, "Table22_Robustness_Proxies.xlsx")
)

# Si tu veux aussi inclure les fichiers EDGE / robustesse supplémentaire,
# décommente les lignes ci-dessous :
# source_files <- c(
#   source_files,
#   file.path(OUT_DIR, "EDGE_Complete_Analysis.xlsx"),
#   file.path(OUT_DIR, "Robustesse_5pct_Results.xlsx")
# )

# Vérifier les fichiers existants
existing_files <- source_files[file.exists(source_files)]

if (length(existing_files) == 0) {
  stop("Aucun fichier Excel source trouvé dans OUT_DIR.")
}

cat("\n[INFO] Fichiers trouvés pour fusion :\n")
print(basename(existing_files))

# Fonction pour noms d’onglets propres et uniques
make_safe_sheet_name <- function(x, used_names = character()) {
  x <- gsub("[\\[\\]\\*\\?/\\\\:]", "_", x)   # caractères interdits Excel
  x <- gsub("\\s+", "_", x)
  x <- substr(x, 1, 31)                       # max 31 caractères Excel
  
  base <- x
  i <- 1
  while (x %in% used_names) {
    suffix <- paste0("_", i)
    x <- substr(base, 1, 31 - nchar(suffix))
    x <- paste0(x, suffix)
    i <- i + 1
  }
  x
}

# Nouveau workbook final
wb_final <- createWorkbook()
used_sheet_names <- character()

for (f in existing_files) {
  
  file_tag <- tools::file_path_sans_ext(basename(f))
  sheets_in_file <- excel_sheets(f)
  
  for (sh in sheets_in_file) {
    
    # lire la feuille
    dat <- tryCatch(
      read_excel(f, sheet = sh),
      error = function(e) NULL
    )
    
    if (is.null(dat)) next
    
    # nom de feuille final
    if (length(sheets_in_file) == 1) {
      proposed_name <- file_tag
    } else {
      proposed_name <- sh
    }
    
    safe_name <- make_safe_sheet_name(proposed_name, used_sheet_names)
    used_sheet_names <- c(used_sheet_names, safe_name)
    
    addWorksheet(wb_final, safe_name)
    writeData(wb_final, sheet = safe_name, x = as.data.frame(dat))
    
    # mise en forme simple
    if (ncol(dat) > 0) {
      setColWidths(wb_final, sheet = safe_name, cols = 1:ncol(dat), widths = "auto")
    }
    freezePane(wb_final, sheet = safe_name, firstRow = TRUE)
  }
}

# Fichier final unique
final_file <- file.path(OUT_DIR, "Tous_les_tableaux_1_a_22.xlsx")
saveWorkbook(wb_final, final_file, overwrite = TRUE)

cat("\n[OK] Fichier final créé :\n")
cat(final_file, "\n")