# Numerical estimation and diagnostics are contained in this document.
# Input matrices are prepared data, not estimates; every render fits anew.
fits <- vector("list", nrow(manifest))
weight_rows <- vector("list", nrow(manifest))
metric_rows <- vector("list", nrow(manifest))
failures <- list()
cached_file <- NULL
for (j in seq_len(nrow(manifest))) {
m <- manifest[j, ]
if (!identical(cached_file, m$matrix_file)) {
source_Y <- as.matrix(read.csv(file.path(params$inputs, m$matrix_file), header=FALSE))
units <- read.csv(file.path(params$inputs, m$units_file))
periods <- read.csv(file.path(params$inputs, m$periods_file))
cached_file <- m$matrix_file
}
ctrl <- which(units$first_treated_column_zero_based == -1)
tr <- which(units$first_treated_column_zero_based == m$cohort_column-1)
Y <- source_Y[c(ctrl,tr), , drop=FALSE]
N0 <- length(ctrl); N1 <- length(tr); T0 <- m$T0; T1 <- m$T1
stopifnot(N0==m$N0, N1==m$N1, ncol(Y)==T0+T1, all(is.finite(Y)), T0>=2, N0>=2)
tic <- proc.time()[[3]]
fit <- tryCatch(synthdid::synthdid_estimate(Y, N0, T0), error=function(e)e)
elapsed <- proc.time()[[3]]-tic
if (inherits(fit,"error")) {
failures[[length(failures)+1]] <- data.frame(m,error=conditionMessage(fit),
donor_noise_sd=sd(as.numeric(apply(Y[seq_len(N0),seq_len(T0),drop=FALSE],1,diff))))
write.csv(do.call(rbind,failures),"failures.csv",row.names=FALSE)
next
}
saveRDS(fit, file.path(params$fit_archive, paste0(m$fit_id,".rds")))
ww <- attr(fit,"weights"); opts <- attr(fit,"opts")
lambda <- as.numeric(ww$lambda); omega <- as.numeric(ww$omega)
stopifnot(length(lambda)==T0, length(omega)==N0, all(is.finite(lambda)),
min(lambda)>=0, abs(sum(lambda)-1)<1e-10, min(omega)>=0, abs(sum(omega)-1)<1e-10)
# Independent four-cell weighted DiD identity, separate from package assembly.
pre <- seq_len(T0); post <- T0+seq_len(T1)
treated <- N0+seq_len(N1)
control_pre <- as.numeric(crossprod(omega,Y[seq_len(N0),pre,drop=FALSE]))
control_post <- sum(omega*rowMeans(Y[seq_len(N0),post,drop=FALSE]))
direct_att <- mean(Y[treated,post,drop=FALSE])-control_post-
sum(lambda*(colMeans(Y[treated,pre,drop=FALSE])-control_pre))
att_check <- abs(direct_att-as.numeric(fit))
stopifnot(att_check < 1e-9*(1+abs(as.numeric(fit))))
# Frank-Wolfe dual gap for the time-weight objective. The package's default
# stopping rule is based on objective improvement, not this certificate.
A <- Y[seq_len(N0),pre,drop=FALSE]
b <- rowMeans(Y[seq_len(N0),post,drop=FALSE])
A <- sweep(A,2,colMeans(A)); b <- b-mean(b)
grad <- as.numeric(2*crossprod(A,A%*%lambda-b)/N0 + 2*opts$zeta.lambda^2*lambda)
fw_gap <- max(0,sum(lambda*grad)-min(grad))
sigma <- sd(as.numeric(apply(Y[seq_len(N0),pre,drop=FALSE],1,diff)))
lvals <- ww$lambda.vals[is.finite(ww$lambda.vals)]
ovals <- ww$omega.vals[is.finite(ww$omega.vals)]
k_recent <- ceiling(.2*T0)
shortest90 <- T0
for (lo in seq_len(T0)) {
ix <- which(cumsum(lambda[lo:T0])>=.9-1e-14)
if (length(ix)) shortest90 <- min(shortest90,ix[1])
}
for (threshold in c(0,1e-8,1e-6,1e-4)) {
active <- lambda>threshold
ix <- which(active)
runs <- sum(active & !c(FALSE,head(active,-1)))
metric_rows[[j]] <- rbind(metric_rows[[j]],data.frame(
fit_id=m$fit_id, threshold=threshold, support=sum(active), support_fraction=mean(active),
support_runs=runs, support_span=if(length(ix))max(ix)-min(ix)+1 else 0))
}
fits[[j]] <- data.frame(m,
estimate=as.numeric(fit), seconds=elapsed,
exact_support=sum(lambda>0), effective_pre_periods=1/sum(lambda^2),
recent_periods=k_recent, mass_last_20pct=sum(tail(lambda,k_recent)),
uniform_mass_last_20pct=k_recent/T0,
weighted_lag=sum(lambda*(T0+1-pre)),
normalized_mean_position=sum(lambda*(pre-1)/(T0-1)),
shortest_90pct_periods=shortest90, shortest_90pct_fraction=shortest90/T0,
lambda_iterations=length(lvals), omega_iterations=length(ovals),
lambda_iteration_cap=length(lvals)>=opts$max.iter,
omega_iteration_cap=length(ovals)>=opts$max.iter,
lambda_final_improvement=if(length(lvals)>=2)tail(lvals,2)[1]-tail(lvals,1) else NA,
default_min_decrease=opts$min.decrease, zeta_lambda=opts$zeta.lambda,
zeta_omega=opts$zeta.omega, lambda_fw_gap=fw_gap,
lambda_fw_gap_noise_scaled=if(sigma>0)fw_gap/sigma^2 else NA,
simplex_error=max(abs(sum(lambda)-1),abs(sum(omega)-1)), att_identity_error=att_check)
weight_rows[[j]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,
specification=m$specification,cohort_column=m$cohort_column,
adoption_source_period=m$adoption_source_period,
pre_period_index=pre,source_period=periods$source_period[pre],
event_period=pre-T0-1,normalized_position=(pre-1)/(T0-1),
lambda=lambda,exact_nonzero=lambda>0)
if(j%%10==0 || j==nrow(manifest)) {
write_json(list(completed=j,total=nrow(manifest),last_fit=m$fit_id),"progress.json",auto_unbox=TRUE)
write.csv(do.call(rbind,fits[seq_len(j)]),"fits.csv",row.names=FALSE)
write.csv(do.call(rbind,weight_rows[seq_len(j)]),"time-weights.csv",row.names=FALSE)
}
}
fits <- do.call(rbind,fits)
weights <- do.call(rbind,weight_rows)
metrics <- do.call(rbind,metric_rows)
if(!length(failures)) write.csv(data.frame(fit_id=character(),error=character()),"failures.csv",row.names=FALSE)
write.csv(fits,"fits.csv",row.names=FALSE)
write.csv(weights,"time-weights.csv",row.names=FALSE)
write.csv(metrics,"support-thresholds.csv",row.names=FALSE)
write_json(list(package="synthdid",version=as.character(packageVersion("synthdid")),
commit=upstream_commit,call="synthdid::synthdid_estimate(Y, N0, T0)",
archive_sha256="9bf048eb12e5972a2290726136acaacbf6f13471072428354eeeed0298607766",
R=R.version.string,session=capture.output(sessionInfo()),
defaults=lapply(formals(synthdid::synthdid_estimate),function(x)paste(deparse(x),collapse=" "))),
"implementation.json",auto_unbox=TRUE,pretty=TRUE)
