---
title: "SDID time weights: second tranche"
subtitle: "Balanced-sample and earlier-cohort rescues · Under development"
date: "2026-10-19"
format:
  html:
    embed-resources: true
    toc: true
    code-fold: true
    html-math-method: katex
execute:
  warning: false
  message: false
params:
  inputs: "../../analyses/sdid-time-weights-tranche2/inputs"
  r_library: "../../.Rlib"
---

## Population and sample construction

This extension covers the **14 papers** classified as possible balanced-sample
or earlier-cohort rescues in the preceding inventory. The frozen roster has
**418 cohort–outcome blocks across 64 prepared specifications**. Each block uses
the source outcome and treatment definition; outcomes are not imputed, and no
new treatment threshold is introduced. The original collection and previously
reported fits remain unchanged.

A deterministic rule is applied before fitting: use complete treated-cohort and
control units over the entire prepared window when that yields at least two
controls and one treated unit. Otherwise, end immediately before the next
adoption cohort, add later adopters while they remain untreated, and retain
complete units in that window. The start of the prepared pre-period is unchanged.
The resulting split is 315 full-window and 103 shortened-window blocks.

These choices change the sample or horizon. Retention, source specification,
window, later-adopter donor counts and earlier preparation steps are exported.
Raw treatment indicators were checked within each selected window: no exposed
donor or treated-unit reversal remains in those windows. This does not make
reversing treatments globally absorbing or remove concerns about anticipation,
spillovers, sample selection, or the causal interpretation of the original design.
Fiscal Capacity retains its existing mixed-IV classification; E-governance
retains its pending classification. These computations do not decide eligibility.

## Estimation target

The target is the **SDID time-weight component**. In the no-covariate author
implementation this is separable from unit-weight fitting. This extension does
not calculate new effect estimates, unit weights, or bootstrap intervals.
The author time solver and its defaults (including sparsification) are used,
with an independent quadratic-program comparison at the same default penalty.
“Optimized” refers to numerical optimization, not penalty selection.

Let A be the donor pre-outcome matrix and b its mean post-outcome vector, both
centered across donors to profile out the free intercept and divided by the
default noise scale. The objective is mean squared prediction error plus
10^-12 times the squared Euclidean norm of the time weights, on the simplex.
A gap of at most 10^-8 and raw feasibility violation of at most 10^-8 are the
numerical objective checks. Such weak regularization can leave weights poorly
determined even when the objective gap is small. Both QP initializations and
rank diagnostics are retained; high-precision follow-up is reported separately.

### Large donor panels

For tall donor matrices, a QR reduction preserves the complete time objective.
If the centered, unscaled matrix D = [A, b] has QR factor R with r rows, use
sqrt(r/N0) R as the input to the **same author solver**, with the intercept already
profiled out. Its mean-square loss, ridge penalty, Frank–Wolfe step, sparsification
and stopping criterion are unchanged in exact arithmetic. No donor is subsampled.
Twenty direct-author steps are compared with reduced steps for every reduced
block, and a full direct-author fit is checked for each preselected example.
Floating-point differences and checks are exported.

```{r}
#| label: setup
#| include: false
.libPaths(c(params$r_library,.libPaths()))
library(synthdid);library(Matrix);library(osqp);library(jsonlite);library(ggplot2);library(knitr)
stopifnot(as.character(packageVersion('synthdid'))=='0.0.9')
manifest <- read.csv(file.path(params$inputs,'manifest.csv'),stringsAsFactors=FALSE)
manifest$representative <- tolower(manifest$representative)=='true'
stopifnot(nrow(manifest)==418,length(unique(manifest$paper_id))==14)
short <- c('aeraej-105'='Right to counsel','aeraej-120'="Women's political voice",'aeraej-159'='Village officials',
 'aeraej-388'='Market integration','aeraej-59'='Procurement auctions','aeraej-6'='Retail pharmacies',
 'aeraej-97'='E-governance rollout','study-02'='Powerful politicians','study-03'='Powerful politicians: reply',
 'study-09'='Vertical integration','study-13'='Ideology and public organizations','study-17'='Trading nontradables',
 'study-18'='Fiscal capacity (mixed-design)','study-22'='Nigerian oil')
```

```{r}
#| label: numerical-methods
#| include: false
author <- function(D,sigma,intercept=TRUE) {
 z <- synthdid:::sc.weight.fw(D,zeta=1e-6*sigma,intercept=intercept,
    min.decrease=1e-5*sigma,max.iter=100)
 z <- synthdid:::sc.weight.fw(D,zeta=1e-6*sigma,intercept=intercept,
    lambda=synthdid:::sparsify_function(z$lambda),min.decrease=1e-5*sigma,max.iter=10000)
 list(x=as.numeric(z$lambda),raw=as.numeric(z$lambda),iterations=sum(is.finite(z$vals)),
      cap=sum(is.finite(z$vals))>=10000,status='author solver returned',repair=0,violation=0)
}
qp <- function(A,b,start) {
 n <- ncol(A);H <- crossprod(A)/nrow(A)+diag(1e-12,n);g <- as.numeric(crossprod(A,b)/nrow(A))
 s <- osqp::osqp(P=2*H,q=-2*g,A=rbind(rep(1,n),Diagonal(n)),l=c(1,rep(0,n)),
   u=c(1,rep(Inf,n)),pars=osqpSettings(eps_abs=1e-11,eps_rel=1e-11,max_iter=1000000L,
   polishing=TRUE,time_limit=30,verbose=FALSE))
 s@WarmStart(x=start,y=rep(0,n+1));z <- s@Solve();raw <- as.numeric(z$x)
 stopifnot(all(is.finite(raw)),sum(pmax(raw,0))>0)
 x <- pmax(raw,0);x <- x/sum(x)
 list(x=x,raw=raw,iterations=z$info$iter,cap=z$info$iter>=1000000L,
 status=z$info$status,repair=sum(abs(x-raw)),violation=max(abs(sum(raw)-1),pmax(0,-min(raw))))
}
objective <- function(x,A,b) mean((as.numeric(A%*%x)-b)^2)+1e-12*sum(x*x)
gap <- function(x,A,b) {
 g <- as.numeric(2*crossprod(A,as.numeric(A%*%x)-b)/nrow(A)+2e-12*x)
 max(0,sum(x*g)-min(g))
}
run_group <- function(ms) {
 con <- file(file.path(params$inputs,ms$matrix_file[1]),'rb')
 source <- matrix(readBin(con,'double',n=ms$source_N[1]*ms$source_T[1],size=8,endian='little'),
   nrow=ms$source_N[1],byrow=TRUE);close(con)
 periods <- read.csv(file.path(params$inputs,ms$periods_file[1]))
 fits <- weights <- failures <- checks <- list()
 for(j in seq_len(nrow(ms))) {
  m <- ms[j,];selection <- read.csv(file.path(params$inputs,m$selection_file))
  ci <- selection$matrix_row_1based[selection$role=='control']
  ti <- selection$matrix_row_1based[selection$role=='treated']
  T0 <- m$T0;T1 <- m$T1;pre <- seq_len(T0);post <- T0+seq_len(T1)
  Y <- source[ci,seq_len(T0+T1),drop=FALSE]
  stopifnot(length(ci)==m$N0,length(ti)==m$N1,all(is.finite(Y)),
   all(is.finite(source[ti,seq_len(T0+T1),drop=FALSE])))
  sigma <- sd(as.numeric(apply(Y[,pre,drop=FALSE],1,diff)))
  if(!is.finite(sigma)||sigma<=0) {
   failures[[length(failures)+1]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,
    method='all',reason='Zero/nonfinite default noise; no substitute weights');next
  }
  D <- cbind(Y[,pre,drop=FALSE],rowMeans(Y[,post,drop=FALSE]));Dc <- sweep(D,2,colMeans(D))
  A <- Dc[,pre,drop=FALSE]/sigma;b <- Dc[,T0+1]/sigma
  reduced <- nrow(D)>4*ncol(D);shortcheck <- fullcheck <- NA_real_;normcheck <- NA_real_
  if(reduced) {
   qrD <- qr(Dc,LAPACK=TRUE);R <- qr.R(qrD);R <- R[,order(qrD$pivot),drop=FALSE]
   DR <- sqrt(nrow(R)/nrow(D))*R
   normcheck <- max(abs(crossprod(DR)/nrow(DR)-crossprod(Dc)/nrow(D)))/max(1,max(abs(crossprod(Dc)/nrow(D))))
   x1 <- synthdid:::sc.weight.fw(D,1e-6*sigma,min.decrease=0,max.iter=20)$lambda
   x2 <- synthdid:::sc.weight.fw(DR,1e-6*sigma,intercept=FALSE,min.decrease=0,max.iter=20)$lambda
   shortcheck <- max(abs(x1-x2));stopifnot(normcheck<1e-10,shortcheck<1e-7)
  } else DR <- D
  results <- list();rankA <- qr(A,tol=1e-10)$rank
  for(method in c('default','qp_uniform','qp_default_start')) {
   tic <- proc.time()[[3]]
   z <- tryCatch(switch(method,default=author(DR,sigma,!reduced),
    qp_uniform=qp(A,b,rep(1/T0,T0)),qp_default_start=qp(A,b,results$default$x)),error=function(e)e)
   elapsed <- proc.time()[[3]]-tic
   if(inherits(z,'error')) {failures[[length(failures)+1]] <- data.frame(fit_id=m$fit_id,
    paper_id=m$paper_id,method=method,reason=conditionMessage(z));next}
   if(method=='default'&&reduced&&m$representative) {
    direct <- author(D,sigma);fullcheck <- max(abs(direct$x-z$x));stopifnot(fullcheck<1e-6)
   }
   x <- z$x;stopifnot(length(x)==T0,all(is.finite(x)),min(x)>=0,abs(sum(x)-1)<1e-9)
   results[[method]] <- z;base <- results$default$x;g <- gap(x,A,b)
   fits[[length(fits)+1]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,paper=unname(short[m$paper_id]),
    specification=m$specification,adoption_source_period=m$adoption_source_period,N0=m$N0,N1=m$N1,
    T0=T0,T1=T1,representative=m$representative,method=method,objective=objective(x,A,b),gap=g,
    accepted=g<=1e-8&&z$violation<=1e-8,iterations=z$iterations,iteration_cap=z$cap,
    solver_status=z$status,raw_violation=z$violation,repair_l1=z$repair,seconds=elapsed,
    l1_from_default=sum(abs(x-base)),max_change_from_default=max(abs(x-base)),
    donor_prediction_rms_change=sqrt(mean((A%*%(x-base))^2)),support_1e6=sum(x>1e-6),
    effective_periods=1/sum(x*x),mass_recent=sum(tail(x,ceiling(.2*T0))),
    mean_position=sum(x*(pre-1)/(T0-1)),noise_scale=sigma,design_rank=rankA,
    rank_deficient_on_simplex=rankA<T0-1,qr_reduction=reduced,
    sample_rule=m$sample_rule,control_retention=m$control_retention,treated_retention=m$treated_retention,
    later_adopter_donors=m$later_adopter_donors,window_end=m$window_end)
   weights[[length(weights)+1]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,method=method,
    source_period=periods$source_period[pre],event_period=pre-T0-1,lambda=x,raw_lambda=z$raw)
  }
  checks[[length(checks)+1]] <- data.frame(fit_id=m$fit_id,qr_reduction=reduced,
    gram_relative_error=normcheck,direct_20_steps_max_difference=shortcheck,
    direct_full_default_max_difference=fullcheck)
 }
 list(fits=do.call(rbind,fits),weights=do.call(rbind,weights),failures=do.call(rbind,failures),checks=do.call(rbind,checks))
}
```

```{r}
#| label: fit-all
#| include: false
groups <- split(manifest,manifest$matrix_file);answers <- list()
for(first in seq(1,length(groups),by=2)) {
 ids <- first:min(first+1,length(groups))
 batch <- parallel::mclapply(groups[ids],run_group,mc.cores=2,mc.preschedule=FALSE)
 if(any(vapply(batch,inherits,logical(1),'try-error')))stop('A worker failed')
 answers <- c(answers,batch)
 for(field in c('fits','weights','failures','checks')) {
  out <- do.call(rbind,lapply(answers,`[[`,field))
  name <- switch(field,weights='time-weights.csv',checks='reduction-checks.csv',paste0(field,'.csv'))
  if(!is.null(out))write.csv(out,name,row.names=FALSE)
 }
 write_json(list(specifications_completed=max(ids),specifications_total=length(groups),
  returned_vectors=nrow(do.call(rbind,lapply(answers,`[[`,'fits')))), 'progress.json',auto_unbox=TRUE)
}
fits <- do.call(rbind,lapply(answers,`[[`,'fits'));weights <- do.call(rbind,lapply(answers,`[[`,'weights'))
if(!file.exists('failures.csv'))write.csv(data.frame(fit_id=character(),paper_id=character(),method=character(),reason=character()),'failures.csv',row.names=FALSE)
a <- subset(fits,method=='default');q <- subset(fits,method=='qp_uniform');b <- subset(fits,method=='qp_default_start')
init <- do.call(rbind,lapply(intersect(q$fit_id,b$fit_id),function(id) {
 x <- subset(weights,fit_id==id & method=='qp_uniform')$lambda
 y <- subset(weights,fit_id==id & method=='qp_default_start')$lambda
 data.frame(fit_id=id,l1_difference=sum(abs(x-y)),max_difference=max(abs(x-y)),
   objective_difference=abs(q$objective[q$fit_id==id]-b$objective[b$fit_id==id]))
}))
write.csv(init,'initialization-sensitivity.csv',row.names=FALSE)
support <- do.call(rbind,lapply(split(weights,list(weights$fit_id,weights$method),drop=TRUE),function(w)
 do.call(rbind,lapply(c(0,1e-8,1e-6,1e-4),function(c) {
  active <- w$lambda>c;ix <- which(active)
  data.frame(fit_id=w$fit_id[1],method=w$method[1],cutoff=c,support=sum(active),
    runs=sum(active & !c(FALSE,head(active,-1))),span=if(length(ix))max(ix)-min(ix)+1 else 0)
 }))))
write.csv(support,'support-thresholds.csv',row.names=FALSE)
write_json(list(author_package='synthdid',version=as.character(packageVersion('synthdid')),
 author_commit='70c1ce3eac58e28c30b67435ca377bb48baa9b8a',osqp=as.character(packageVersion('osqp')),
 eta_lambda=1e-6,normalized_ridge=1e-12,default_max_iterations=10000,
 default_sparsification=TRUE,objective_gap_target=1e-8,time_weights_only=TRUE,
 penalty_retuning=FALSE,original_results_modified=FALSE), 'implementation.json',pretty=TRUE,auto_unbox=TRUE,digits=16)
```

## Numerical results

```{r}
#| echo: false
kable(aggregate(cbind(accepted,iteration_cap)~method,data=fits,sum))
kable(head(init[order(-init$l1_difference),],12),digits=8)
```

The CSVs retain all zeros and small weights. Support thresholds are diagnostic;
no thresholded vector is substituted for a fitted vector. Multiple outcomes and
cohorts share data and donor pools, so the blocks are dependent examples.

- [Sample manifest](manifest.csv)
- [Fit diagnostics](fits.csv) / [all vectors](time-weights.csv)
- [Initializer sensitivity](initialization-sensitivity.csv)
- [QR equivalence checks](reduction-checks.csv)
- [Input provenance](input-provenance.json)
- [Executable source](analysis.qmd)