R/RcppExports.R

In coresynth: Fast and Unified Synthetic Control Methods

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' Fast Interactive Fixed Effects (IFE) for Generalized Synthetic Control
#'
#' Implements Xu (2017) IFE model with optional covariate adjustment.
#' When X_co has p > 0 slices, runs an EM loop alternating between:
#'   E-step: truncated SVD of Y_tilde = Y_co - X_co * beta
#'   M-step: panel OLS to update beta given current factors
#' When X_co has 0 slices (default), falls back to the plain 3-step estimator.
#'
#' @param Y_co     Control units outcome matrix (T x N_co)
#' @param Y_tr_pre Treated units pre-treatment outcomes (T_pre x N_tr)
#' @param r        Number of latent factors (must be <= min(T, N_co))
#' @param X_co     Time-varying covariate cube (T x N_co x p). Pass an empty
#'                 cube (0 slices) for the covariate-free estimator.
#' @param X_tr_pre Time-varying covariate cube for treated units in the
#'                 pre-treatment window (T_pre x N_tr x p). Required for
#'                 correct Step 2 loading estimation per Xu (2017): lambda_hat
#'                 is estimated from Y_tr_pre - X_tr_pre * beta (covariate-
#'                 demeaned). Pass an empty cube (0 slices) to skip demeaning
#'                 (backward-compatible, but biased when beta != 0).
#' @param max_iter Maximum EM iterations (default 50)
#' @param tol      Convergence tolerance on relative beta change (default 1e-6)
#' @return A list with components:
#'   * `F`: estimated time factors (T x r).
#'   * `L_co`: control-unit factor loadings (N_co x r).
#'   * `L_tr`: treated-unit factor loadings (N_tr x r).
#'   * `Y_tr_hat`: estimated treated-unit counterfactual outcomes (T x N_tr).
#'   * `singular_values`: singular values from the final truncated SVD.
#'   * `beta`: estimated covariate coefficients (p x 1), empty when no
#'     covariates are supplied.
#' @export
gsc_ife_cpp <- function(Y_co, Y_tr_pre, r, X_co, X_tr_pre, max_iter = 50L, tol = 1e-6) {
    .Call(`_coresynth_gsc_ife_cpp`, Y_co, Y_tr_pre, r, X_co, X_tr_pre, max_iter, tol)
}

#' Fast Placebo Test for SDID
#'
#' For each control unit, treats it as the "pseudo-treated" unit and
#' estimates the leave-one-out SDID effect. The distribution of these
#' placebo effects provides a permutation-based null distribution for inference.
#'
#' @param Y_pre     Control units pre-treatment outcomes (T_pre x N_co)
#' @param Y_post    Control units post-treatment outcomes (T_post x N_co)
#' @param time_weights Lambda weights for pre-treatment periods (T_pre x 1)
#' @param zeta2     Ridge penalty (same as used in the main estimate)
#' @return A numeric vector of length `N_co`. Each element is the
#'   leave-one-out placebo SDID effect obtained by treating that control unit
#'   as the pseudo-treated unit; the vector serves as a permutation-based null
#'   distribution for inference.
#' @export
sdid_placebo_cpp <- function(Y_pre, Y_post, time_weights, zeta2) {
    .Call(`_coresynth_sdid_placebo_cpp`, Y_pre, Y_post, time_weights, zeta2)
}

#' Fast Leave-One-Out Placebo Test for SCM (Abadie et al. 2010)
#'
#' For each control unit, treats it as pseudo-treated and fits SCM weights
#' from the remaining N_co-1 donors. Returns MSPE components for constructing
#' MSPE-ratio permutation p-values in R.
#'
#' @param Y_pre   Control pre-treatment outcomes (T_pre x N_co)
#' @param Y_post  Control post-treatment outcomes (T_post x N_co)
#' @param max_iter Outer coordinate-descent iterations (default 100)
#' @param tol      Convergence tolerance for V updates (default 1e-4)
#' @param z_rows  Optional 1-based pre-period row indices of the outer
#'   evaluation window (the `v_window` of the treated fit), so each placebo
#'   refit optimises V on the same window. `NULL` (default) uses all rows.
#'   MSPE components are always computed on the full pre/post windows.
#' @param wolfe If `TRUE`, each placebo refit uses the Wolfe min-norm-point
#'   inner solver, matching a treated fit made with `qp_solver = "wolfe"` so
#'   the permutation stays symmetric.
#' @return A list with:
#'   * `mspe_pre`:  N_co-vector of pre-treatment MSPE per placebo unit
#'   * `mspe_post`: N_co-vector of post-treatment MSPE per placebo unit
#'   * `effects`:   N_co-vector of mean post-period gap per placebo unit
#'   * `gaps`:      (T_pre + T_post) x N_co matrix of placebo gap paths
#' @export
scm_placebo_cpp <- function(Y_pre, Y_post, max_iter = 100L, tol = 1e-4, z_rows = NULL, wolfe = FALSE) {
    .Call(`_coresynth_scm_placebo_cpp`, Y_pre, Y_post, max_iter, tol, z_rows, wolfe)
}

#' Fast Leave-One-Out Placebo Test for SCM with a Predictor Specification
#'
#' Covariate-spec counterpart of [scm_placebo_cpp()]: for each control unit,
#' treats it as pseudo-treated with its own predictor column `X0[, i]` and
#' fits the nested V/W optimisation against the remaining donors' predictors
#' `X0[, -i]`, evaluating the prediction loss on pre-treatment outcomes.
#' Each leave-one-out problem is identical to a [scm_weights_cpp()] call on
#' the same submatrices; iterations are independent and run in parallel
#' under OpenMP.
#'
#' @param X0      Predictor matrix for control units (k x N_co), on the same
#'   scale as the treated fit (SD-scaled when `scale_predictors = TRUE`)
#' @param Y_pre   Control pre-treatment outcomes (T_pre x N_co)
#' @param Y_post  Control post-treatment outcomes (T_post x N_co)
#' @param max_iter Outer coordinate-descent iterations (default 100)
#' @param tol      Convergence tolerance for V updates (default 1e-4)
#' @param z_rows  Optional 1-based pre-period row indices of the outer
#'   evaluation window (the `v_window` of the treated fit), so each placebo
#'   refit optimises V on the same window. `NULL` (default) uses all rows.
#'   MSPE components are always computed on the full pre/post windows.
#' @param multistart If `TRUE`, each placebo refit uses the same
#'   deterministic multi-start outer search as the treated fit, keeping the
#'   permutation test symmetric.
#' @param wolfe If `TRUE`, each placebo refit uses the Wolfe min-norm-point
#'   inner solver, matching a treated fit made with `qp_solver = "wolfe"`.
#' @return A list with:
#'   * `mspe_pre`:  N_co-vector of pre-treatment MSPE per placebo unit
#'   * `mspe_post`: N_co-vector of post-treatment MSPE per placebo unit
#'   * `effects`:   N_co-vector of mean post-period gap per placebo unit
#'   * `gaps`:      (T_pre + T_post) x N_co matrix of placebo gap paths
#'   A placebo unit whose solver fails yields NaN entries.
#' @export
scm_placebo_x_cpp <- function(X0, Y_pre, Y_post, max_iter = 100L, tol = 1e-4, z_rows = NULL, multistart = FALSE, wolfe = FALSE) {
    .Call(`_coresynth_scm_placebo_x_cpp`, X0, Y_pre, Y_post, max_iter, tol, z_rows, multistart, wolfe)
}

#' Fast Matrix Completion using Soft-Impute Algorithm
#'
#' Solves: min_L (1/2) ||O o (Y - L)||_F^2 + lambda * ||L||_*
#' via iterative SVD soft-thresholding (Mazumder, Hastie, Tibshirani 2010).
#' Note: lambda is NOT normalized by |O|. Default lambda = 0.01 * sigma_max(Y).
#'
#' @param Y       Observed outcome matrix (N x T). Unobserved entries should be 0.
#' @param O       Binary mask matrix (N x T): 1 = observed, 0 = missing (treated post).
#' @param lambda  Nuclear norm penalty (soft-threshold on singular values).
#' @param max_iter Maximum iterations.
#' @param tol     Convergence tolerance (relative Frobenius norm change).
#' @return A numeric matrix of the same dimension as `Y` (N x T): the
#'   completed low-rank matrix `L` that minimises the soft-thresholded
#'   nuclear-norm objective.
#' @export
soft_impute_cpp <- function(Y, O, lambda, max_iter = 1000L, tol = 1e-5) {
    .Call(`_coresynth_soft_impute_cpp`, Y, O, lambda, max_iter, tol)
}

proj_simplex <- function(y) {
    .Call(`_coresynth_proj_simplex`, y)
}

solve_simplex_qp <- function(Q, c, max_iter = 10000L, tol = 1e-6, x0 = NULL) {
    .Call(`_coresynth_solve_simplex_qp`, Q, c, max_iter, tol, x0)
}

#' SCM Inner Weights (QP Given V)
#'
#' Solves the inner-loop QP for SCM: given a fixed diagonal metric matrix V,
#' finds donor weights W on the simplex minimising the V-weighted covariate
#' loss. The returned weights are a KKT-verified exact optimum whenever the
#' active-set solver converges, with accelerated projected gradient as a
#' fallback.
#'
#' @param X0     Covariate matrix for control units (k x N_co)
#' @param X1     Covariate vector for the treated unit (k x 1)
#' @param V_diag Diagonal of the metric matrix V (k x 1, non-negative, need not sum to 1)
#' @param wolfe  If `TRUE`, solve with the Wolfe min-norm-point method, which
#'   returns a Caratheodory-sparse optimum (at most k+1 donors carry weight)
#'   instead of one arbitrary point of a degenerate optimal face. `FALSE`
#'   (default) uses the warm-started active-set solver.
#' @return Donor weight vector W (N_co x 1) on the unit simplex
#' @export
scm_inner_weights_cpp <- function(X0, X1, V_diag, wolfe = FALSE) {
    .Call(`_coresynth_scm_inner_weights_cpp`, X0, X1, V_diag, wolfe)
}

#' SCM Outer Weights (Joint Optimization of W and V)
#'
#' Jointly optimises donor weights W (on the simplex) and the diagonal
#' metric matrix V via coordinate descent on the pre-treatment prediction
#' MSPE, following Abadie, Diamond & Hainmueller (2010).
#'
#' When `t_train > 0`, V is selected by minimising MSPE on a validation
#' window (rows t_train..T_pre-1 of Z) while W is fitted on the full
#' predictor matrix X. This is appropriate when X is a fixed predictor
#' matrix that contains no validation-period outcome information (the
#' user-supplied predictors case). For the outcomes-only case the proper
#' Abadie (2021) S.3.2 train/validation split is implemented in R
#' (`.scm_oos_outcomes()`): candidate W(V) are fitted on training-half
#' outcomes only, by passing the training rows as X and the validation
#' rows as Z to this function with `t_train = -1`.
#'
#' @param X0      Covariate matrix for control units (k x N_co, typically pre-treatment outcomes)
#' @param X1      Covariate vector for the treated unit (k x 1)
#' @param Z0      Outcome matrix for control units in the pre-treatment window (T_pre x N_co)
#' @param Z1      Outcome vector for the treated unit in the pre-treatment window (T_pre x 1)
#' @param max_iter Maximum coordinate-descent iterations (default 100)
#' @param tol     Convergence tolerance on MSPE improvement (default 1e-4)
#' @param t_train Validation-window split for V selection.
#'   -1 (default): V selected on the full Z window (in-sample).
#'   Positive: rows t_train..(T_pre-1) of Z form the validation window used
#'   to select V (W is fitted on the full X throughout); after selecting V*,
#'   W is refit and the reported loss uses the full Z window.
#' @param z_rows Optional 1-based row indices of Z defining the evaluation
#'   window for the outer V optimisation (the `v_window` argument of
#'   [scm_fit()]). `NULL` (default) evaluates on the full Z window. Takes
#'   precedence over `t_train`; the reported loss always uses the full Z
#'   window.
#' @param multistart If `TRUE`, the outer V optimisation runs a
#'   deterministic multi-start search (screened start set, coordinate-descent
#'   polish, Nelder-Mead refinement) instead of a single coordinate-descent
#'   pass from the uniform V. The result is never worse (in outer loss) than
#'   the single-start path.
#' @param cheap_face If `TRUE`, the inner simplex QP first tries a cheap
#'   bordered-KKT direct solve on each active-set face, falling back to the
#'   scale-robust null-space solve when that system is singular. Worth trying
#'   only in the outcomes-only regime (no user predictors, no predictor
#'   rescaling), where V is dense and the faces that carry the fit are well
#'   conditioned; it reproduces the null-space solution to round-off.
#'   `FALSE` (default) goes straight to the null-space solve, which is what
#'   predictor fits need for scale invariance and rank-deficient faces.
#' @param wolfe If `TRUE`, every inner QP is solved with the Wolfe
#'   min-norm-point method, which returns a Caratheodory-sparse optimum (at
#'   most k+1 donors carry weight) instead of one arbitrary point of a
#'   degenerate optimal face. `FALSE` (default) uses the warm-started
#'   active-set solver.
#' @return A list with:
#'   * `W`: Donor weight vector (N_co x 1) on the unit simplex
#'   * `V`: Optimal metric diagonal (k x 1, normalised to sum to 1)
#'   * `loss`: Final pre-treatment prediction loss (full pre-treatment window)
#' @export
scm_weights_cpp <- function(X0, X1, Z0, Z1, max_iter = 100L, tol = 1e-4, t_train = -1L, z_rows = NULL, multistart = FALSE, cheap_face = FALSE, wolfe = FALSE) {
    .Call(`_coresynth_scm_weights_cpp`, X0, X1, Z0, Z1, max_iter, tol, t_train, z_rows, multistart, cheap_face, wolfe)
}

#' Calculate SDID Unit Weights (omega)
#'
#' Solves the regularized QP:
#' min over omega in Delta: sum_t (sum_i omega_i Y_it - Y_tr_t)^2 + zeta^2 * T_pre * ||omega||^2
#'
#' This corresponds to equation (5) in Arkhangelsky et al. (2021).
#'
#' @param Y_pre Pre-treatment outcome matrix for control units (T_pre x N_co)
#' @param Y_tr_pre Pre-treatment outcome vector for treated unit (T_pre x 1), averaged if multiple
#' @param zeta2 Ridge penalty parameter (zeta^2). The code internally multiplies by T_pre per the paper.
#' @return A numeric vector of length `N_co` holding the SDID unit weights
#'   `omega` (non-negative and summing to one).
#' @export
sdid_unit_weights_cpp <- function(Y_pre, Y_tr_pre, zeta2) {
    .Call(`_coresynth_sdid_unit_weights_cpp`, Y_pre, Y_tr_pre, zeta2)
}

#' Calculate SDID Time Weights (lambda)
#'
#' Solves the time-weight QP (with implicit intercept lambda_0 concentrated out):
#'
#' min over lambda in Delta_pre: ||Y_post_target - Y_pre_co^T lambda||^2 + zeta_t^2 * N_co * ||lambda||^2
#'
#' The caller is responsible for pre-demeaning Y_pre_co (row-wise) and
#' Y_post_target (subtract the cross-unit mean) to concentrate out lambda_0,
#' as described in Arkhangelsky et al. (2021) Algorithm 1, Eq. (2.3).
#'
#' @param Y_pre_co  Pre-treatment outcomes for control units, row-demeaned (T_pre x N_co)
#' @param Y_post_target Post-treatment mean per control unit, demeaned (N_co x 1)
#' @param zeta_t    Ridge penalty for time weights (paper: 1e-6 * sigma_hat)
#' @return A numeric vector of length `T_pre` holding the SDID time weights
#'   `lambda` (non-negative and summing to one).
#' @export
sdid_time_weights_cpp <- function(Y_pre_co, Y_post_target, zeta_t) {
    .Call(`_coresynth_sdid_time_weights_cpp`, Y_pre_co, Y_post_target, zeta_t)
}

#' Calculate SDID Estimate (tau_sdid)
#'
#' Given unit weights omega and time weights lambda, computes the SDID
#' estimator as a weighted two-way difference:
#'
#' tau_sdid = (Y_tr_post_mean - Y_tr_pre_wt) - (Y_co_post_wt - Y_co_pre_wt)
#'
#' @param Y_pre_co   Control pre-treatment outcomes (T_pre x N_co)
#' @param Y_post_co  Control post-treatment outcomes (T_post x N_co)
#' @param Y_pre_tr   Treated pre-treatment outcomes (T_pre x 1)
#' @param Y_post_tr  Treated post-treatment outcomes (T_post x 1)
#' @param omega      Unit weights (N_co x 1)
#' @param lambda     Time weights (T_pre x 1)
#' @return A single numeric value: the SDID treatment-effect estimate
#'   `tau_sdid`.
#' @export
sdid_estimate_cpp <- function(Y_pre_co, Y_post_co, Y_pre_tr, Y_post_tr, omega, lambda) {
    .Call(`_coresynth_sdid_estimate_cpp`, Y_pre_co, Y_post_co, Y_pre_tr, Y_post_tr, omega, lambda)
}

#' Tensor Unfolding (Matricization) for Synthetic Interventions
#'
#' @param T_cube A 3D array (cube) of dimensions (n1, n2, n3)
#' @param mode The mode to unfold along (1, 2, or 3)
#' @return A numeric matrix: the mode-`mode` unfolding (matricization) of
#'   `T_cube`, with dimensions `n1 x (n2 * n3)`, `n2 x (n1 * n3)`, or
#'   `n3 x (n1 * n2)` for `mode` 1, 2, or 3 respectively.
#' @export
tensor_unfold_cpp <- function(T_cube, mode) {
    .Call(`_coresynth_tensor_unfold_cpp`, T_cube, mode)
}

#' SI-PCR: Synthetic Interventions via Principal Component Regression
#'
#' Implements the SI-PCR estimator of Agarwal et al. (2025).
#' Uses the top-k SVD of pre-treatment control outcomes to find donor
#' weights that predict each treated unit's pre-treatment trajectory,
#' then applies those weights to post-treatment control outcomes.
#'
#' @param Y_pre_co  Pre-treatment control outcomes (T_pre x N_co)
#' @param Y_post_co Post-treatment control outcomes (T_post x N_co)
#' @param Y_pre_tr  Pre-treatment treated outcomes (T_pre x N_tr)
#' @param k         Number of SVD components to retain
#' @return A list with:
#'   * `W`: Donor weight matrix (N_co x N_tr)
#'   * `Y_hat`: Counterfactual post-treatment outcomes (T_post x N_tr)
#' @export
si_pcr_cpp <- function(Y_pre_co, Y_post_co, Y_pre_tr, k) {
    .Call(`_coresynth_si_pcr_cpp`, Y_pre_co, Y_post_co, Y_pre_tr, k)
}

#' Kalman Filter and RTS Smoother (TASC)
#'
#' Implements the Kalman filter (forward pass) and Rauch-Tung-Striebel smoother
#' (backward pass) for the state-space model in Rho et al. (2026):
#'
#'   State:       z(t+1) = A z(t) + C + eta(t),  eta(t) ~ N(0, Q)
#'   Observation: y_t    = W z_t + eps_t,        eps_t  ~ N(0, R)
#'
#' Observation rows with NA (treated post-intervention) are automatically
#' dropped at each time step so only control-unit rows update the filter.
#'
#' The P update uses the numerically stable Joseph form:
#'   P(t|t) = (I - K W_obs) P(t|t-1) (I - K W_obs)^T + K R_obs K^T
#'
#' @param Y  Observed data matrix (N x T). Use NA for unobserved entries.
#' @param W  Observation / loading matrix (N x r)
#' @param A  State transition matrix (r x r). Pass diag(r) for random-walk dynamics.
#' @param C  State drift vector (r x 1)
#' @param Q  State noise covariance (r x r)
#' @param R  Observation noise covariance (N x N, diagonal in practice)
#' @param z0 Initial state mean (r x 1)
#' @param P0 Initial state covariance (r x r)
#' @return A list with z_smooth, P_smooth, P_cross, z_pred, z_upd.
#'   P_cross is an r x r x (T-1) cube. Slice t (C++ 0-indexed, t=0,...,T-2)
#'   stores P(t+1, t | T) (0-indexed), i.e. P(t+2, t+1 | T) in 1-indexed
#'   Shumway-Stoffer notation. Formula: P(t+1|T) * J_t^T (eq. 6.68-6.69).
#' @export
kalman_smoother_cpp <- function(Y, W, A, C, Q, R, z0, P0) {
    .Call(`_coresynth_kalman_smoother_cpp`, Y, W, A, C, Q, R, z0, P0)
}