R/erri.R

In ERRI: Economic Resilience and Recovery Index

#' Estimate the Economic Resilience and Recovery Index
#'
#' Constructs a counterfactual outcome path from pre-shock observations and
#' measures the magnitude and duration of the deviation after the shock.
#'
#' @param data A data frame containing a regularly ordered time variable and
#'   numeric outcome.
#' @param time Character string naming the time column.
#' @param outcome Character string naming the outcome column.
#' @param shock_time Shock date or time. For grouped data, either one common
#'   value or a named vector with names equal to the group values.
#' @param unit Optional character string naming a grouping column.
#' @param method Counterfactual method: linear time `"trend"`, pre-shock
#'   `"mean"`, or autoregressive `"ar1"`.
#' @param scale Scale for gaps: pre-shock standard deviation (`"sd"`), absolute
#'   pre-shock mean (`"mean"`), or no scaling (`"none"`).
#' @param epsilon Recovery tolerance in scaled-gap units.
#' @param consecutive Number of consecutive observations within `epsilon`
#'   required to declare recovery.
#' @param weights Named non-negative weights for `resistance`, `loss`,
#'   `recovery`, `strength`, `stability`, and `transformation`.
#' @param level Confidence level for counterfactual prediction intervals.
#'
#' @return An object of class `erri` containing component estimates, scores,
#'   the composite index, trajectories, settings, and fitted models.
#' @export
#' @examples
#' dat <- erri_example_data()
#' fit <- erri(dat, time = "year", outcome = "income",
#'             unit = "region", shock_time = 2020)
#' fit
erri <- function(data, time, outcome, shock_time, unit = NULL,
                 method = c("trend", "mean", "ar1"),
                 scale = c("sd", "mean", "none"), epsilon = 0.25,
                 consecutive = 2L, weights = NULL, level = 0.95) {
    method <- match.arg(method)
    scale <- match.arg(scale)
    if (!is.data.frame(data)) .erri_stop("'data' must be a data frame.")
    needed <- c(time, outcome, unit)
    if (!all(needed %in% names(data))) {
        .erri_stop("These columns were not found: ",
                   paste(setdiff(needed, names(data)), collapse = ", "), ".")
    }
    if (!is.numeric(data[[outcome]])) .erri_stop("The outcome must be numeric.")
    if (!is.numeric(epsilon) || length(epsilon) != 1L || epsilon < 0) {
        .erri_stop("'epsilon' must be one non-negative number.")
    }
    consecutive <- as.integer(consecutive)
    if (length(consecutive) != 1L || is.na(consecutive) || consecutive < 1L) {
        .erri_stop("'consecutive' must be a positive integer.")
    }
    if (!is.numeric(level) || length(level) != 1L || level <= 0 || level >= 1) {
        .erri_stop("'level' must be between zero and one.")
    }
    weights <- .validate_weights(weights)
    groups <- if (is.null(unit)) rep("All", nrow(data)) else as.character(data[[unit]])
    ids <- unique(groups)
    components <- vector("list", length(ids))
    trajectories <- vector("list", length(ids))
    models <- vector("list", length(ids))
    names(trajectories) <- names(models) <- ids

    for (i in seq_along(ids)) {
        id <- ids[i]
        d <- data[groups == id, , drop = FALSE]
        d <- d[order(.as_time_numeric(d[[time]])), , drop = FALSE]
        if (anyDuplicated(d[[time]])) .erri_stop("Duplicated times in unit '", id, "'.")
        st <- if (length(shock_time) == 1L) shock_time else {
            if (is.null(names(shock_time)) || !id %in% names(shock_time)) {
                .erri_stop("Grouped 'shock_time' must be scalar or named by unit.")
            }
            shock_time[[id]]
        }
        sidx <- .match_shock(d[[time]], st)
        if (sidx <= 5L || sidx > nrow(d)) {
            .erri_stop("Shock for unit '", id,
                       "' must follow at least five pre-shock observations.")
        }
        y <- d[[outcome]]
        cf <- .fit_counterfactual(d[[time]], y, sidx, method, level)
        sc <- .scale_value(y[seq_len(sidx - 1L)], scale)
        gap <- (cf$fit - y) / sc
        post <- sidx:length(y)
        pgap <- gap[post]
        adverse <- pmax(pgap, 0)
        depth <- if (all(is.na(adverse))) NA_real_ else max(adverse, na.rm = TRUE)
        loss <- if (all(is.na(adverse))) NA_real_ else sum(adverse, na.rm = TRUE)
        rtime <- .run_recovery(pgap, epsilon, consecutive)
        trough <- if (all(is.na(adverse))) 1L else which.max(adverse)
        end_idx <- if (is.na(rtime)) length(pgap) else max(rtime + 1L, trough)
        strength <- if (end_idx <= trough || !is.finite(depth)) 0 else {
            pmax(depth - adverse[end_idx], 0) / (end_idx - trough)
        }
        pre_resid <- y[seq_len(sidx - 1L)] - cf$fit[seq_len(sidx - 1L)]
        post_resid <- y[post] - cf$fit[post]
        denom <- stats::sd(pre_resid, na.rm = TRUE)
        stability <- stats::sd(post_resid, na.rm = TRUE) /
            ifelse(is.finite(denom) && denom > 0, denom, sc)
        trans_start <- if (is.na(rtime)) length(post) else min(rtime + 1L, length(post))
        transformation <- mean(pmax(-pgap[trans_start:length(pgap)], 0), na.rm = TRUE)
        if (!is.finite(transformation)) transformation <- 0
        components[[i]] <- data.frame(
            unit = id, shock_time = as.character(d[[time]][sidx]),
            shock_depth = depth, cumulative_loss = loss,
            recovery_time = if (is.na(rtime)) length(post) else rtime,
            recovered = !is.na(rtime), recovery_strength = strength,
            stability_ratio = stability, transformation = transformation,
            stringsAsFactors = FALSE
        )
        trajectories[[id]] <- data.frame(
            unit = id, time = d[[time]], observed = y,
            counterfactual = cf$fit, lower = cf$lower, upper = cf$upper,
            scaled_gap = gap, post_shock = seq_along(y) >= sidx,
            stringsAsFactors = FALSE
        )
        models[[id]] <- cf$model
    }
    comp <- do.call(rbind, components)
    scores <- .component_scores(comp)
    comp <- cbind(comp, setNames(100 * scores, paste0(names(scores), "_score")))
    comp$ERRI <- 100 * as.numeric(as.matrix(scores) %*% weights)
    rownames(comp) <- NULL
    out <- list(results = comp, trajectories = trajectories, models = models,
                weights = weights,
                settings = list(time = time, outcome = outcome, unit = unit,
                                method = method, scale = scale,
                                epsilon = epsilon, consecutive = consecutive,
                                level = level, shock_time = shock_time),
                call = match.call())
    class(out) <- "erri"
    out
}

#' @export
print.erri <- function(x, digits = 2, ...) {
    cat("Economic Resilience and Recovery Index (ERRI)\n")
    cat("Counterfactual:", x$settings$method, "| Scale:", x$settings$scale, "\n\n")
    keep <- c("unit", "shock_depth", "cumulative_loss", "recovery_time",
              "recovered", "stability_ratio", "transformation", "ERRI")
    print(x$results[keep], digits = digits, row.names = FALSE)
    invisible(x)
}

#' @export
summary.erri <- function(object, ...) {
    ans <- list(call = object$call, settings = object$settings,
                weights = object$weights, results = object$results)
    class(ans) <- "summary.erri"
    ans
}