R/cross_sectional.R

In EventStudy: Event Study Analysis

#' Cross-Sectional Regression of CARs
#'
#' Regress cumulative abnormal returns (CARs) on firm characteristics to
#' explain cross-sectional variation in event effects. Supports OLS with
#' heteroskedasticity-consistent (HC) standard errors.
#'
#' @param task A fitted EventStudyTask with abnormal returns computed.
#' @param formula A formula with the response on the left (ignored; CAR is
#'   always the dependent variable) and explanatory variables on the right,
#'   e.g., \code{~ log_market_cap + leverage}.
#' @param data A data frame of firm characteristics. Must contain an
#'   \code{event_id} column to merge with CARs.
#' @param car_window Two-element integer vector specifying the CAR window
#'   as \code{c(start, end)} relative indices. Default is the full event window.
#' @param robust Logical. If TRUE and the \pkg{sandwich} package is available,
#'   compute HC1 robust standard errors. Default TRUE.
#'
#' @return A list with class \code{"es_cross_sectional"} containing:
#'   \describe{
#'     \item{model}{The fitted \code{lm} object}
#'     \item{coefficients}{Coefficient table with (robust) standard errors}
#'     \item{r_squared}{R-squared of the regression}
#'     \item{n_obs}{Number of observations}
#'     \item{car_data}{The merged CAR + characteristics data}
#'   }
#'
#' @family eventstudy-statistics
#' @export
cross_sectional_regression <- function(task, formula, data,
                                        car_window = NULL,
                                        robust = TRUE) {
  if (!inherits(task, "EventStudyTask")) {
    rlang::abort(
      paste0("`task` must be an EventStudyTask object, not <", class(task)[1], ">."),
      class = c("eventstudy_error_bad_argument", "eventstudy_error"))
  }
  if (!"model" %in% names(task$data_tbl)) {
    rlang::abort(
      "Abnormal returns not computed. Run fit_model() first.",
      class = c("eventstudy_error_not_fitted", "eventstudy_error"))
  }
  if (!"event_id" %in% names(data)) {
    rlang::abort(
      "`data` must contain an `event_id` column for merging.",
      class = c("eventstudy_error_missing_column", "eventstudy_error"))
  }

  # A5 (2026-09-24): duplicated event_id in `data` silently corrupts the
  # join (many-to-many inflation) -- error instead of proceeding.
  dup_ids <- unique(data$event_id[duplicated(data$event_id)])
  if (length(dup_ids) > 0) {
    rlang::abort(
      paste0("`data` contains duplicated event_id value(s): ",
             paste(dup_ids, collapse = ", ")),
      class = c("eventstudy_error_bad_argument", "eventstudy_error"))
  }

  # A5: car_window extending outside the available event window of ANY
  # event is an error (not a silently-empty or NA CAR).
  if (!is.null(car_window)) {
    event_ranges <- purrr::map2_dfr(
      task$data_tbl$event_id, task$data_tbl$data,
      function(eid, d) {
        ew <- d %>% dplyr::filter(event_window == 1)
        tibble::tibble(
          event_id = eid,
          min_ri = suppressWarnings(min(ew$relative_index, na.rm = TRUE)),
          max_ri = suppressWarnings(max(ew$relative_index, na.rm = TRUE))
        )
      }
    )
    out_of_range <- event_ranges %>%
      dplyr::filter(car_window[1] < min_ri | car_window[2] > max_ri)
    if (nrow(out_of_range) > 0) {
      finite_ranges <- event_ranges %>%
        dplyr::filter(is.finite(min_ri), is.finite(max_ri))
      avail_min <- suppressWarnings(min(finite_ranges$min_ri))
      avail_max <- suppressWarnings(max(finite_ranges$max_ri))
      rlang::abort(
        paste0("car_window = c(", car_window[1], ", ", car_window[2],
               ") extends outside the available event window [",
               avail_min, ", ", avail_max, "] for event_id(s): ",
               paste(out_of_range$event_id, collapse = ", ")),
        class = c("eventstudy_error_bad_argument", "eventstudy_error"))
    }
  }

  # Extract CARs for each event
  cars <- .extract_cars(task, car_window)

  # A5: task events without a matching row in `data` are dropped with one
  # warning (before the inner_join silently drops them).
  unmatched <- setdiff(cars$event_id, data$event_id)
  if (length(unmatched) > 0) {
    warning(
      "cross_sectional_regression: ", length(unmatched),
      " task event(s) have no matching row in `data` and were dropped ",
      "(event_id: ", paste(unmatched, collapse = ", "), ").",
      call. = FALSE
    )
  }

  # Merge with firm characteristics
  merged <- dplyr::inner_join(cars, data, by = "event_id")

  if (nrow(merged) == 0) {
    stop("No matching event_id values between task and data.")
  }

  # A5: an event whose CAR is NA (a missing AR inside car_window --
  # .extract_cars() no longer na.rm's the window sum) is excluded from the
  # regression with one warning listing the event_id(s), rather than being
  # silently included as an NA row (which lm() would drop unexplained).
  na_car_events <- merged$event_id[is.na(merged$car)]
  if (length(na_car_events) > 0) {
    warning(
      "cross_sectional_regression: ", length(na_car_events),
      " event(s) excluded because an abnormal return inside car_window is ",
      "missing (event_id: ", paste(na_car_events, collapse = ", "), ").",
      call. = FALSE
    )
    merged <- merged %>% dplyr::filter(!is.na(car))
  }

  if (nrow(merged) == 0) {
    stop("No matching event_id values between task and data.")
  }

  # Build regression formula: car ~ <user formula RHS>
  rhs <- as.character(formula)[length(as.character(formula))]
  reg_formula <- stats::as.formula(paste("car ~", rhs))

  # Fit OLS -- wrap in tryCatch so an lm() error (e.g., singular formula)
  # returns NULL + warning instead of crashing the caller.
  fit <- tryCatch(
    stats::lm(reg_formula, data = merged),
    error = function(e) {
      warning("cross_sectional_regression: lm() failed -- ",
              conditionMessage(e), call. = FALSE)
      NULL
    }
  )
  if (is.null(fit)) return(NULL)
  fit_summary <- summary(fit)

  # Rank-deficiency produces NA coefficients for aliased columns; vcovHC on
  # a rank-deficient matrix calls solve() which throws a singular error.  Warn
  # the caller before attempting the sandwich path so the message is always
  # emitted regardless of whether sandwich is installed.
  rank_deficient <- fit$rank < length(stats::coef(fit))
  if (rank_deficient) {
    warning(
      "cross_sectional_regression: design matrix is rank-deficient; ",
      "aliased columns dropped by lm(). ",
      "Robust SEs from vcovHC may be unreliable on near-singular designs.",
      call. = FALSE
    )
  }

  # Robust standard errors
  if (robust && requireNamespace("sandwich", quietly = TRUE)) {
    # Wrap vcovHC -- on a perfectly singular design solve() throws here.
    vcov_hc <- tryCatch(
      sandwich::vcovHC(fit, type = "HC1"),
      error = function(e) {
        warning("cross_sectional_regression: vcovHC() failed (singular design) -- ",
                "returning NA standard errors.", call. = FALSE)
        NULL
      }
    )
    all_coefs <- stats::coef(fit)   # length = all terms including aliased (NA) ones
    all_names <- names(all_coefs)
    if (!is.null(vcov_hc)) {
      # vcovHC drops aliased (NA) coefficients, so its row/col names are a
      # subset of all_names.  Build a full-length SE vector with NA in the
      # aliased positions so the data.frame rows align correctly.
      se_fitted <- sqrt(diag(vcov_hc))
      se <- setNames(rep(NA_real_, length(all_names)), all_names)
      se[names(se_fitted)] <- se_fitted
    } else {
      se <- rep(NA_real_, length(all_names))
    }
    coef_tbl <- data.frame(
      estimate  = all_coefs,
      std.error = se,
      statistic = all_coefs / se,
      p.value   = 2 * stats::pt(abs(all_coefs / se),
                                  df = fit$df.residual, lower.tail = FALSE),
      row.names = all_names
    )
  } else {
    coef_tbl <- as.data.frame(fit_summary$coefficients)
    names(coef_tbl) <- c("estimate", "std.error", "statistic", "p.value")
    if (robust) {
      # Upgrade from message() to warning() and name the lost capability and
      # the OLS fallback that is being used instead (EXTERNAL-03 pipeline half).
      warning(
        "Package 'sandwich' not available: cluster/robust SEs unavailable. ",
        "Falling back to OLS standard errors. ",
        "Install 'sandwich' for heteroskedasticity-consistent SEs.",
        call. = FALSE
      )
    }
  }

  result <- list(
    model = fit,
    coefficients = coef_tbl,
    r_squared = fit_summary$r.squared,
    adj_r_squared = fit_summary$adj.r.squared,
    n_obs = nrow(merged),
    car_data = merged
  )
  class(result) <- "es_cross_sectional"
  result
}


#' Extract CARs from Task
#' @noRd
.extract_cars <- function(task, car_window) {
  purrr::pmap_dfr(
    list(task$data_tbl$event_id, task$data_tbl$firm_symbol,
         task$data_tbl$group, task$data_tbl$data),
    function(eid, sym, grp, d) {
      event_data <- d %>%
        dplyr::filter(event_window == 1)

      if (!"abnormal_returns" %in% names(event_data)) {
        return(tibble::tibble(
          event_id = eid, firm_symbol = sym, group = grp, car = NA_real_
        ))
      }

      if (!is.null(car_window)) {
        event_data <- event_data %>%
          dplyr::filter(relative_index >= car_window[1],
                        relative_index <= car_window[2])
      }

      # A5 (2026-09-24): NO na.rm -- ANY missing AR inside the window makes
      # the CAR for that event NA (not the sum of the finite subset). This
      # is the shared .extract_cars() convention inherited by
      # cross_sectional_regression() (which excludes+warns), car_by_group(),
      # car_quantiles() and plot_car_distribution() (which silently propagate
      # the NA per their existing na.rm=TRUE summary statistics).
      ar_vals <- event_data$abnormal_returns
      car_val <- if (length(ar_vals) == 0) NA_real_ else sum(ar_vals)
      tibble::tibble(
        event_id = eid,
        firm_symbol = sym,
        group = grp,
        car = car_val
      )
    }
  )
}


#' Format Cross-Sectional Regression Results
#'
#' Builds the character vector rendered by \code{print.es_cross_sectional}
#' (one element per output line). The coefficient table is rendered via
#' \code{base::print.data.frame} and spliced in at the same position.
#'
#' @param x An \code{es_cross_sectional} object.
#' @param ... Additional arguments (unused).
#'
#' @return A character vector, one element per printed line.
#'
#' @export
format.es_cross_sectional <- function(x, ...) {
  utils::capture.output({
    cat("Cross-Sectional Regression of CARs\n")
    cat("===================================\n")
    cat("N:", x$n_obs, "\n")
    cat("R-squared:", round(x$r_squared, 4), "\n")
    cat("Adj. R-squared:", round(x$adj_r_squared, 4), "\n\n")
    cat("Coefficients:\n")
    base::print.data.frame(round(x$coefficients, 6))
  })
}

#' Print Cross-Sectional Regression Results
#'
#' @param x An \code{es_cross_sectional} object.
#' @param ... Additional arguments (unused).
#'
#' @export
print.es_cross_sectional <- function(x, ...) {
  cat(format(x), sep = "\n")
  invisible(x)
}


#' Compare CARs Across Groups
#'
#' Performs a t-test (two groups) or ANOVA (3+ groups) to test whether
#' CARs differ significantly across groups.
#'
#' @param task A fitted EventStudyTask.
#' @param group_var Name of the grouping variable. Defaults to "group"
#'   (the standard group column in EventStudyTask).
#' @param car_window Optional two-element vector for CAR window.
#'
#' @return A list with test results and group-level summary statistics.
#'
#' @family eventstudy-statistics
#' @export
car_by_group <- function(task, group_var = "group", car_window = NULL) {
  cars <- .extract_cars(task, car_window)

  if (!group_var %in% names(cars)) {
    rlang::abort(
      paste0("`group_var` = \"", group_var, "\" not found."),
      class = c("eventstudy_error_unknown_column", "eventstudy_error"))
  }

  groups <- unique(cars[[group_var]])
  n_groups <- length(groups)

  # Summary by group
  group_summary <- cars %>%
    dplyr::group_by(.data[[group_var]]) %>%
    dplyr::summarise(
      n = dplyr::n(),
      mean_car = mean(car, na.rm = TRUE),
      sd_car = sd(car, na.rm = TRUE),
      median_car = stats::median(car, na.rm = TRUE),
      min_car = min(car, na.rm = TRUE),
      max_car = max(car, na.rm = TRUE),
      .groups = "drop"
    )

  # Test
  if (n_groups == 2) {
    g1 <- cars$car[cars[[group_var]] == groups[1]]
    g2 <- cars$car[cars[[group_var]] == groups[2]]
    test <- stats::t.test(g1, g2)
    test_name <- "Welch Two-Sample t-test"
  } else if (n_groups > 2) {
    fml <- stats::as.formula(paste("car ~", group_var))
    test <- stats::oneway.test(fml, data = cars)
    test_name <- "One-way ANOVA (Welch)"
  } else {
    test <- NULL
    test_name <- "Only one group -- no test performed"
  }

  list(
    summary = group_summary,
    test = test,
    test_name = test_name,
    n_groups = n_groups
  )
}


#' CAR Quantiles
#'
#' Compute quantiles of the CAR distribution across events.
#'
#' @param task A fitted EventStudyTask.
#' @param probs Probability vector for quantiles. Default
#'   \code{c(0.05, 0.25, 0.5, 0.75, 0.95)}.
#' @param car_window Optional two-element vector for CAR window.
#'
#' @return A named numeric vector of quantiles.
#'
#' @family eventstudy-statistics
#' @export
car_quantiles <- function(task, probs = c(0.05, 0.25, 0.5, 0.75, 0.95),
                           car_window = NULL) {
  cars <- .extract_cars(task, car_window)
  stats::quantile(cars$car, probs = probs, na.rm = TRUE)
}


#' Plot CAR Distribution
#'
#' Create a histogram of CARs across events, optionally colored by group.
#'
#' @param task A fitted EventStudyTask.
#' @param car_window Optional two-element vector for CAR window.
#' @param bins Number of histogram bins. Default 30.
#' @param by_group Logical. If TRUE, color histogram by group.
#' @param title Optional plot title.
#'
#' @return A ggplot2 object.
#'
#' @family eventstudy-plots
#' @export
plot_car_distribution <- function(task, car_window = NULL,
                                   bins = 30, by_group = FALSE,
                                   title = NULL) {
  cars <- .extract_cars(task, car_window)

  if (is.null(title)) title <- "Distribution of Cumulative Abnormal Returns"

  p <- ggplot2::ggplot(cars, ggplot2::aes(x = car))

  if (by_group && length(unique(cars$group)) > 1) {
    p <- p +
      ggplot2::geom_histogram(ggplot2::aes(fill = group),
                               bins = bins, alpha = 0.7,
                               position = "identity") +
      ggplot2::scale_fill_brewer(palette = "Set2")
  } else {
    p <- p +
      ggplot2::geom_histogram(fill = "steelblue", alpha = 0.7, bins = bins)
  }

  p +
    ggplot2::geom_vline(xintercept = 0, linetype = "dashed", color = "red") +
    ggplot2::labs(title = title, x = "CAR", y = "Count") +
    ggplot2::theme_minimal() +
    ggplot2::theme(plot.title = ggplot2::element_text(hjust = 0.5))
}