R/counterfactual.R

In ERRI: Economic Resilience and Recovery Index

.fit_counterfactual <- function(time, y, shock_idx,
                                method = c("trend", "mean", "ar1"),
                                level = 0.95) {
    method <- match.arg(method)
    pre <- seq_len(shock_idx - 1L)
    if (length(pre) < 5L) {
        .erri_stop("At least five pre-shock observations are required.")
    }
    tt <- .as_time_numeric(time)
    tt <- (tt - tt[1L]) / max(diff(range(tt)), 1)
    alpha <- 1 - level

    if (method == "mean") {
        fit <- mean(y[pre], na.rm = TRUE)
        pred <- rep(fit, length(y))
        res <- y[pre] - fit
        se <- stats::sd(res, na.rm = TRUE) * sqrt(1 + 1 / sum(is.finite(y[pre])))
        lower <- pred + stats::qnorm(alpha / 2) * se
        upper <- pred + stats::qnorm(1 - alpha / 2) * se
        model <- list(method = "mean", fitted = rep(fit, length(pre)),
                      residuals = res, pre = pre)
    } else if (method == "trend") {
        d <- data.frame(y = y[pre], time = tt[pre])
        model_fit <- stats::lm(y ~ time, data = d, na.action = stats::na.exclude)
        pp <- stats::predict(model_fit, newdata = data.frame(time = tt),
                             interval = "prediction", level = level)
        pred <- pp[, "fit"]
        lower <- pp[, "lwr"]
        upper <- pp[, "upr"]
        model <- list(method = "trend", fit = model_fit,
                      fitted = stats::fitted(model_fit),
                      residuals = stats::residuals(model_fit), pre = pre)
    } else {
        yy <- y[pre]
        lag_y <- yy[-length(yy)]
        now_y <- yy[-1L]
        model_fit <- stats::lm(now_y ~ lag_y)
        cf <- rep(NA_real_, length(y))
        cf[pre] <- c(yy[1L], stats::fitted(model_fit))
        phi <- stats::coef(model_fit)[2L]
        intercept <- stats::coef(model_fit)[1L]
        last <- yy[length(yy)]
        for (j in shock_idx:length(y)) {
            last <- intercept + phi * last
            cf[j] <- last
        }
        sig <- stats::sd(stats::residuals(model_fit), na.rm = TRUE)
        h <- pmax(seq_along(y) - shock_idx + 1L, 0L)
        # A direct geometric sum remains non-negative for stationary,
        # unit-root, and explosive AR(1) estimates. The common closed form
        # can suffer a negative denominator and produce NaN when abs(phi) > 1.
        se <- sig * sqrt(vapply(pmax(h, 1L), function(k) {
            sum(phi^(2 * seq.int(0L, k - 1L)))
        }, numeric(1)))
        pred <- cf
        lower <- pred + stats::qnorm(alpha / 2) * se
        upper <- pred + stats::qnorm(1 - alpha / 2) * se
        model <- list(method = "ar1", fit = model_fit,
                      fitted = cf[pre], residuals = stats::residuals(model_fit),
                      pre = pre)
    }
    list(fit = as.numeric(pred), lower = as.numeric(lower),
         upper = as.numeric(upper), model = model)
}