diff --git a/R/cv.R b/R/cv.R index f9ad35e..3bceddc 100644 --- a/R/cv.R +++ b/R/cv.R @@ -308,6 +308,33 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix message("Some units have too few pre-treatment observations. Remove them automatically in Cross-Validation.\n") } + ## ----------------------------------------------------## + ## Pre-compute fold-specific training data (speed) + ## ----------------------------------------------------## + ## These objects depend only on the CV split, not on r/lambda. + ## Precomputing avoids repeated TT*N matrix copies inside inner loops. + II.cv.list <- vector("list", k) + YY.cv.list <- vector("list", k) + W.cv.list <- NULL + if (use_weight == 1) { + W.cv.list <- vector("list", k) + } + for (ii in 1:k) { + II.cv.tmp <- II + II.cv.tmp[rmCV[[ii]]] <- 0 + II.cv.list[[ii]] <- II.cv.tmp + + YY.cv.tmp <- YY + YY.cv.tmp[rmCV[[ii]]] <- 0 + YY.cv.list[[ii]] <- YY.cv.tmp + + if (use_weight == 1) { + W.tmp <- W.use + W.tmp[rmCV[[ii]]] <- 0 + W.cv.list[[ii]] <- W.tmp + } + } + ## get count matrix if (use_weight == 0) { count.T.cv <- count.T.cv.old <- table(T.on) @@ -359,6 +386,17 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix CV.out.ife[, "r"] <- c(r.old:r.max) CV.out.ife[, "PC"] <- CV.out.ife[, "GMoment"] <- CV.out.ife[, "Moment"] <- CV.out.ife[, "MAD"] <- CV.out.ife[, "MSPE"] <- CV.out.ife[, "WMSPE"] <- CV.out.ife[, "GMSPE"] <- CV.out.ife[, "WGMSPE"] <- 1e20 + ## warm-start buffers (per fold) to reduce iterations across r + fit.start.ife <- vector("list", k) + beta.start.ife <- vector("list", k) + for (ii in 1:k) { + fit.start.ife[[ii]] <- as.matrix(Y0CV[, , ii]) + beta.start.ife[[ii]] <- as.matrix(beta0CV[, , ii]) + } + ## warm-start buffer for the full-sample fit across r + fit.start.ife.full <- as.matrix(Y0) + beta.start.ife.full <- as.matrix(beta0) + for (i in 1:dim(CV.out.ife)[1]) { ## cross-validation loop starts ## inter FE based on control, before & after r <- CV.out.ife[i, "r"] @@ -373,21 +411,28 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix index.moment.list <- c() MAD.list <- c() for (ii in 1:k) { - II.cv <- II - II.cv[rmCV[[ii]]] <- 0 - YY.cv <- YY - YY.cv[rmCV[[ii]]] <- 0 if (use_weight) { - W.use2 <- W.use - W.use2[rmCV[[ii]]] <- 0 + W.use2 <- W.cv.list[[ii]] } else { W.use2 <- as.matrix(0) } - est.cv.fit <- inter_fe_ub( - YY.cv, as.matrix(Y0CV[, , ii]), X, II.cv, - W.use2, as.matrix(beta0CV[, , ii]), + est.cv.obj <- inter_fe_ub( + YY.cv.list[[ii]], + fit.start.ife[[ii]], + X, + II.cv.list[[ii]], + W.use2, + beta.start.ife[[ii]], r, force, tol, max.iteration - )$fit + ) + est.cv.fit <- est.cv.obj$fit + ## update warm-starts for next r on the same fold + fit.start.ife[[ii]] <- est.cv.fit + if (p > 0 && !is.null(est.cv.obj$beta)) { + beta.tmp <- est.cv.obj$beta + beta.tmp[is.na(beta.tmp)] <- 0 + beta.start.ife[[ii]] <- as.matrix(beta.tmp) + } index.cv <- as.character(T.on[estCV[[ii]]]) index.cv[which(is.na(index.cv))] <- "Control" weight.cv <- count.T.cv[index.cv] @@ -453,18 +498,25 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix gmoment <- sum(weight.cv.g * resid.g.mean) / sum(weight.cv) } + ## overall fit (warm-start across r) est.cv <- inter_fe_ub( YY, - Y0, + fit.start.ife.full, X, II, W.use, - beta0, + beta.start.ife.full, r, force, tol, max.iteration - ) ## overall + ) + fit.start.ife.full <- est.cv$fit + if (p > 0 && !is.null(est.cv$beta)) { + beta.tmp <- est.cv$beta + beta.tmp[is.na(beta.tmp)] <- 0 + beta.start.ife.full <- as.matrix(beta.tmp) + } sigma2 <- est.cv$sigma2 IC <- est.cv$IC PC <- est.cv$PC @@ -666,6 +718,17 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix break_count <- 0 break_check <- 0 + ## warm-start buffers (per fold) to reduce iterations across lambda + fit.start.mc <- vector("list", k) + beta.start.mc <- vector("list", k) + for (ii in 1:k) { + fit.start.mc[[ii]] <- as.matrix(Y0CV[, , ii]) + beta.start.mc[[ii]] <- as.matrix(beta0CV[, , ii]) + } + ## warm-start buffer for the full-sample fit across lambda + fit.start.mc.full <- as.matrix(Y0) + beta.start.mc.full <- as.matrix(beta0) + for (i in 1:length(lambda)) { ## k <- 5 SSE <- 0 @@ -678,21 +741,28 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix MAD.list <- c() for (ii in 1:k) { - II.cv <- II - II.cv[rmCV[[ii]]] <- 0 - YY.cv <- YY - YY.cv[rmCV[[ii]]] <- 0 if (use_weight) { - W.use2 <- W.use - W.use2[rmCV[[ii]]] <- 0 + W.use2 <- W.cv.list[[ii]] } else { W.use2 <- as.matrix(0) } - est.cv.fit <- inter_fe_mc( - YY.cv, as.matrix(Y0CV[, , ii]), - X, II.cv, W.use2, as.matrix(beta0CV[, , ii]), + est.cv.obj <- inter_fe_mc( + YY.cv.list[[ii]], + fit.start.mc[[ii]], + X, + II.cv.list[[ii]], + W.use2, + beta.start.mc[[ii]], 1, lambda[i], force, tol, max.iteration - )$fit + ) + est.cv.fit <- est.cv.obj$fit + ## update warm-starts for next lambda on the same fold + fit.start.mc[[ii]] <- est.cv.fit + if (p > 0 && !is.null(est.cv.obj$beta)) { + beta.tmp <- est.cv.obj$beta + beta.tmp[is.na(beta.tmp)] <- 0 + beta.start.mc[[ii]] <- as.matrix(beta.tmp) + } index.cv <- as.character(T.on[estCV[[ii]]]) index.cv[which(is.na(index.cv))] <- "Control" weight.cv <- count.T.cv[index.cv] @@ -750,11 +820,23 @@ fect_cv <- function(Y, # Outcome variable, (T*N) matrix moment <- sum(weight.cv * resid.mean) / sum(weight.cv) gmoment <- sum(weight.cv.g * resid.g.mean) / sum(weight.cv) + ## overall fit (warm-start across lambda) est.cv <- inter_fe_mc( - YY, Y0, X, II, W.use, beta0, + YY, + fit.start.mc.full, + X, + II, + W.use, + beta.start.mc.full, 1, lambda[i], force, tol, max.iteration - ) ## overall + ) + fit.start.mc.full <- est.cv$fit + if (p > 0 && !is.null(est.cv$beta)) { + beta.tmp <- est.cv$beta + beta.tmp[is.na(beta.tmp)] <- 0 + beta.start.mc.full <- as.matrix(beta.tmp) + } eff.v.cv <- c(Y - est.cv$fit)[cv.pos] meff <- as.numeric(tapply(eff.v.cv, t.on.cv, mean))