R/shocks.R

In ERRI: Economic Resilience and Recovery Index

#' Screen a series for a single structural shock
#'
#' Evaluates admissible split points using the proportional reduction in
#' residual sum of squares from separate linear trends. This is a screening
#' diagnostic, not a causal identification procedure.
#'
#' @param data A data frame.
#' @param time Character string naming the ordered time column.
#' @param outcome Character string naming the numeric outcome column.
#' @param min_segment Minimum observations on each side of a candidate split.
#' @param direction `"any"`, `"down"`, or `"up"`.
#' @return A data frame of candidate times, scores, and estimated level shifts,
#'   sorted from strongest to weakest.
#' @export
#' @examples
#' dat <- subset(erri_example_data(), region == "North")
#' head(detect_shocks(dat, "year", "income"))
detect_shocks <- function(data, time, outcome, min_segment = 5L,
                          direction = c("any", "down", "up")) {
    direction <- match.arg(direction)
    if (!is.data.frame(data) || !all(c(time, outcome) %in% names(data))) {
        .erri_stop("Supply a data frame and valid 'time' and 'outcome' columns.")
    }
    if (!is.numeric(data[[outcome]])) .erri_stop("The outcome must be numeric.")
    min_segment <- as.integer(min_segment)
    d <- data[order(.as_time_numeric(data[[time]])), , drop = FALSE]
    n <- nrow(d)
    if (is.na(min_segment) || min_segment < 3L || n < 2L * min_segment) {
        .erri_stop("Not enough observations for the requested segment length.")
    }
    tt <- .as_time_numeric(d[[time]])
    yy <- d[[outcome]]
    full <- stats::lm(yy ~ tt)
    rss0 <- sum(stats::residuals(full)^2, na.rm = TRUE)
    candidates <- min_segment:(n - min_segment)
    ans <- lapply(candidates, function(k) {
        left_dat <- data.frame(y = yy[seq_len(k)], x = tt[seq_len(k)])
        right_idx <- (k + 1L):n
        right_dat <- data.frame(y = yy[right_idx], x = tt[right_idx])
        left <- stats::lm(y ~ x, data = left_dat)
        right <- stats::lm(y ~ x, data = right_dat)
        rss <- sum(stats::residuals(left)^2, na.rm = TRUE) +
            sum(stats::residuals(right)^2, na.rm = TRUE)
        before <- stats::predict(left, newdata = data.frame(x = tt[k]))
        after <- stats::predict(right, newdata = data.frame(x = tt[k + 1L]))
        shift <- unname(after - before)
        score <- pmax((rss0 - rss) / pmax(rss0, .Machine$double.eps), 0)
        data.frame(candidate_time = as.character(d[[time]][k + 1L]),
                   score = score, level_shift = shift)
    })
    out <- do.call(rbind, ans)
    if (direction == "down") out <- out[out$level_shift < 0, , drop = FALSE]
    if (direction == "up") out <- out[out$level_shift > 0, , drop = FALSE]
    out <- out[order(out$score, decreasing = TRUE), , drop = FALSE]
    rownames(out) <- NULL
    out
}