---
title: "Synthetic DiD time weights: defaults versus tighter optimization"
subtitle: "Fixed tuning parameters and samples · exploratory methodological check"
date: "2026-10-19"
format:
  html:
    embed-resources: true
    toc: true
    code-fold: true
execute:
  warning: false
  message: false
params:
  inputs: "../../analyses/sdid-time-weights-default/inputs"
  r_library: "../../.Rlib"
---

## What is being compared

This exercise holds the outcome matrices, donor pools, cohort windows,
intercept and **default time-weight penalty** fixed. “Optimized” means numerical
optimization of that same objective, **not cross-validation of a new penalty**.
It concerns time weights only: unit weights, ATTs and inference are not refit or
promoted to the main collection. In the default no-covariate model the two weight
problems are separate, so improving unit-weight convergence is not necessary to
study the time-weight solution.

Four algorithm settings are compared for every available block:

1. **Author default:** reproduce the authors' time-weight solver, including its
   100-iteration preliminary stage, sparsification, 10,000-iteration final cap
   and `min.decrease = 1e-5 * noise.level`.
2. **Author tighter, sparsification:** same algorithm, but a 100,000 final cap and
   `min.decrease = 1e-8 * noise.level`. The preliminary cap remains 100.
3. **Author tighter, no sparsification:** a single uniform-start author-solver
   run with that tighter stopping rule and 100,000 cap.
4. **QP optimized:** an independent full-simplex quadratic program for the
   identical objective, solved with OSQP. Uniform and default-weight initializers
   are both retained to expose numerical sensitivity in nearly flat problems.

The pinned author package is `synth-inference/synthdid` 0.0.9, commit
`70c1ce3eac58e28c30b67435ca377bb48baa9b8a`. The default penalty coefficient is
`zeta.lambda = 1e-6 * noise.level`. After dividing the objective by the positive
noise variance, the ridge coefficient is therefore **1e-12**. Very similar
objective values need not imply similar weights at such weak regularization.

```{r}
#| label: setup
#| include: false
.libPaths(c(params$r_library,.libPaths()))
library(synthdid);library(Matrix);library(osqp);library(ggplot2);library(jsonlite);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'
old <- read.csv('../sdid-time-weights-default/time-weights.csv')
short <- read.csv('../sdid-time-weights-default/paper-summary.csv')
short_names <- setNames(short$paper,short$paper_id)
methods <- c('default','author_tight_sparse','author_tight_no_sparse','qp_uniform','qp_default_start')
method_labels <- c(default='Author default',author_tight_sparse='Author tighter + sparsify',
  author_tight_no_sparse='Author tighter, no sparsify',qp_uniform='QP optimized',
  qp_default_start='QP: default initializer')
gap_tolerance <- 1e-8
constraint_tolerance <- 1e-8
support_cutoffs <- c(0,1e-8,1e-6,1e-4)
```

## Objective and acceptance checks

Let A contain the controls' pre-period outcomes, with each column centered across
controls, and b their centered mean post-period outcome. The program minimizes
`mean((A %*% lambda - b)^2) + zeta.lambda^2 * sum(lambda^2)` subject to nonnegative
weights summing to one. Centering profiles out the authors' free intercept.

For comparability across outcomes the objective and gradient are divided by
`noise.level^2`. A simplex Frank–Wolfe gap, `lambda' gradient - min(gradient)`,
is an upper bound on objective suboptimality in exact arithmetic. The numerical
acceptance target is **gap ≤ 1e-8**, with simplex feasibility and no material
negative weights. Solver status alone is not the acceptance criterion.

OSQP uses absolute and relative tolerances 1e-11, a 1,000,000 iteration limit,
polishing, and a 30-second limit per solve. Its raw weights and feasibility
violations are preserved. For feasible weight comparisons only, negative
roundoff is clipped to zero and the vector is normalized; the repair magnitude
and the post-repair gap are checked. Tiny positive QP numbers are not claimed
to be mathematically nonzero. The main support comparison uses **weight > 1e-6**,
with three other cutoffs provided.

```{r}
#| label: numerical-functions
#| include: false
objective <- function(x,A,b) mean((as.numeric(A%*%x)-b)^2)+1e-12*sum(x*x)
dual_gap <- function(x,A,b) {
  grad <- as.numeric(2*crossprod(A,as.numeric(A%*%x)-b)/nrow(A)+2e-12*x)
  max(0,sum(x*grad)-min(grad))
}
author_fit <- function(D,sigma,tighter=FALSE,sparse=TRUE) {
  tol <- (if(tighter)1e-8 else 1e-5)*sigma
  limit <- if(tighter)100000L else 10000L
  a <- synthdid:::sc.weight.fw(D,zeta=1e-6*sigma,intercept=TRUE,
    min.decrease=tol,max.iter=if(sparse)100L else limit)
  if(sparse) a <- synthdid:::sc.weight.fw(D,zeta=1e-6*sigma,intercept=TRUE,
    lambda=synthdid:::sparsify_function(a$lambda),min.decrease=tol,max.iter=limit)
  list(x=as.numeric(a$lambda),raw=as.numeric(a$lambda),iterations=sum(is.finite(a$vals)),
    limit=limit,status='returned',repair_l1=0,raw_violation=0)
}
qp_fit <- function(A,b,start) {
  n <- ncol(A);H <- crossprod(A)/nrow(A)+diag(1e-12,n)
  g <- as.numeric(crossprod(A,b)/nrow(A))
  solver <- 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))
  solver@WarmStart(x=start,y=rep(0,n+1))
  z <- solver@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,limit=1000000L,status=z$info$status,
    repair_l1=sum(abs(x-raw)),raw_violation=max(abs(sum(raw)-1),max(0,-min(raw))))
}
fit_group <- function(rows) {
  m0 <- rows[1,]
  source <- as.matrix(read.csv(file.path(params$inputs,m0$matrix_file),header=FALSE))
  units <- read.csv(file.path(params$inputs,m0$units_file))
  periods <- read.csv(file.path(params$inputs,m0$periods_file))
  control <- source[units$first_treated_column_zero_based==-1,,drop=FALSE]
  out <- list();fail <- list();vectors <- list();support <- list()
  for(j in seq_len(nrow(rows))) {
    m <- rows[j,];T0 <- m$T0;pre <- seq_len(T0);post <- T0+seq_len(m$T1)
    D <- cbind(control[,pre,drop=FALSE],rowMeans(control[,post,drop=FALSE]))
    sigma <- sd(as.numeric(apply(control[,pre,drop=FALSE],1,diff)))
    if(!is.finite(sigma)||sigma==0) {
      fail[[length(fail)+1]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,
        reason='Zero/nonfinite default noise scale; no substitute time weights assigned')
      next
    }
    A <- sweep(control[,pre,drop=FALSE],2,colMeans(control[,pre,drop=FALSE]))/sigma
    b <- (rowMeans(control[,post,drop=FALSE])-mean(control[,post,drop=FALSE]))/sigma
    results <- list()
    for(method in methods) {
      tic <- proc.time()[[3]]
      z <- tryCatch(switch(method,
        default=author_fit(D,sigma),
        author_tight_sparse=author_fit(D,sigma,TRUE,TRUE),
        author_tight_no_sparse=author_fit(D,sigma,TRUE,FALSE),
        qp_uniform=qp_fit(A,b,rep(1/T0,T0)),
        qp_default_start=qp_fit(A,b,results$default$x)),error=function(e)e)
      seconds <- proc.time()[[3]]-tic
      if(inherits(z,'error')) {
        fail[[length(fail)+1]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,
          reason=paste(method,conditionMessage(z)));next
      }
      x <- z$x;stopifnot(length(x)==T0,all(is.finite(x)),min(x)>=-1e-12,abs(sum(x)-1)<1e-9)
      if(method=='default') {
        prior <- old$lambda[old$fit_id==m$fit_id]
        stopifnot(length(prior)==T0,max(abs(x-prior))<1e-11)
      }
      results[[method]] <- z
      base <- results$default$x;gap <- dual_gap(x,A,b)
      out[[length(out)+1]] <- data.frame(fit_id=m$fit_id,paper_id=m$paper_id,
        paper=unname(short_names[m$paper_id]),specification=m$specification,
        adoption_source_period=m$adoption_source_period,N0=m$N0,N1=m$N1,T0=T0,T1=m$T1,
        representative=m$representative,method=method,objective=objective(x,A,b),
        gap=gap,accepted=gap<=gap_tolerance && z$raw_violation<=constraint_tolerance,
        iterations=z$iterations,iteration_cap=z$iterations>=z$limit,solver_status=z$status,
        raw_violation=z$raw_violation,repair_l1=z$repair_l1,seconds=seconds,
        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)))
      vectors[[length(vectors)+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)
      for(cutoff in support_cutoffs) {
        active <- x>cutoff;ix <- which(active)
        support[[length(support)+1]] <- data.frame(fit_id=m$fit_id,method=method,
          cutoff=cutoff,support=sum(active),runs=sum(active & !c(FALSE,head(active,-1))),
          span=if(length(ix))max(ix)-min(ix)+1 else 0)
      }
    }
  }
  list(fits=do.call(rbind,out),weights=do.call(rbind,vectors),
    support=do.call(rbind,support),failures=do.call(rbind,fail))
}
```

```{r}
#| label: fit-all
#| include: false
# Only prepared input matrices are loaded. Author defaults are recomputed;
# old CSVs are used only to validate equality with the previously delivered run.
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],fit_group,mc.cores=2,mc.preschedule=FALSE)
  if(any(vapply(batch,inherits,logical(1),'try-error')))stop('A worker failed')
  answers <- c(answers,batch)
  write_json(list(specifications_completed=max(ids),specifications_total=length(groups)),
    'progress.json',auto_unbox=TRUE)
}
fits <- do.call(rbind,lapply(answers,`[[`,'fits'))
weights <- do.call(rbind,lapply(answers,`[[`,'weights'))
support <- do.call(rbind,lapply(answers,`[[`,'support'))
failures <- do.call(rbind,lapply(answers,`[[`,'failures'))
write.csv(fits,'fits.csv',row.names=FALSE)
write.csv(weights,'time-weights.csv',row.names=FALSE)
write.csv(support,'support-thresholds.csv',row.names=FALSE)
write.csv(failures,'failures.csv',row.names=FALSE)
write_json(list(author_package=as.character(packageVersion('synthdid')),
  author_commit='70c1ce3eac58e28c30b67435ca377bb48baa9b8a',osqp=as.character(packageVersion('osqp')),
  time_penalty_multiplier=1e-6,gap_tolerance=gap_tolerance,constraint_tolerance=constraint_tolerance,
  author_tight_max_iter=100000,author_tight_min_decrease_multiplier=1e-8,
  samples_changed=FALSE,penalties_retuned=FALSE,unit_weights_refit=FALSE,
  session=capture.output(sessionInfo())),'implementation.json',auto_unbox=TRUE,pretty=TRUE)
```

## Numerical results

```{r}
#| label: summaries
#| echo: false
method_summary <- do.call(rbind,lapply(methods,function(method){
  f <- fits[fits$method==method,]
  data.frame(method=method,n=nrow(f),accepted=sum(f$accepted),iteration_caps=sum(f$iteration_cap),
    median_gap=median(f$gap),max_gap=max(f$gap),median_l1_change=median(f$l1_from_default),
    max_l1_change=max(f$l1_from_default))
}))
write.csv(method_summary,'method-summary.csv',row.names=FALSE)
kable(method_summary,digits=8)
def <- fits[fits$method=='default',];q <- fits[fits$method=='qp_uniform',]
alt <- fits[fits$method=='qp_default_start',]
comparison <- merge(def,q,by='fit_id',suffixes=c('_default','_qp'))
comparison$objective_reduction <- comparison$objective_default-comparison$objective_qp
comparison$support_change <- comparison$support_1e6_qp-comparison$support_1e6_default
comparison$recent_mass_change <- comparison$mass_recent_qp-comparison$mass_recent_default
write.csv(comparison,'default-vs-optimized.csv',row.names=FALSE)
initialization <- do.call(rbind,lapply(q$fit_id,function(id){
  a <- weights$lambda[weights$fit_id==id&weights$method=='qp_uniform']
  b <- weights$lambda[weights$fit_id==id&weights$method=='qp_default_start']
  data.frame(fit_id=id,l1_difference=sum(abs(a-b)),max_difference=max(abs(a-b)),
    both_accepted=q$accepted[q$fit_id==id]&&alt$accepted[alt$fit_id==id],
    objective_difference=abs(q$objective[q$fit_id==id]-alt$objective[alt$fit_id==id]))
}))
write.csv(initialization,'initialization-sensitivity.csv',row.names=FALSE)
```

The acceptance count measures a numerical objective-gap criterion, not statistical
validity or identification of the weights. With the default tiny ridge, substantial
weight changes can occur along nearly flat directions. Initializer comparisons
and changes in donor predictions are therefore reported alongside objective gains.

## Representative comparisons

The representative cohort/outcome for each paper is the same one selected before
the earlier default run; no example is chosen for having a large difference.

```{r}
#| label: profiles
#| echo: false
#| fig-width: 12
#| fig-height: 13
representatives <- def[def$representative,]
w <- weights[weights$fit_id%in%representatives$fit_id & weights$method%in%c('default','qp_uniform'),]
w$method <- factor(w$method,levels=c('default','qp_uniform'),labels=c('Author default','QP optimized'))
w$paper <- unname(short_names[w$paper_id])
p <- ggplot(w,aes(event_period,lambda,color=method,shape=method))+
  geom_point(size=1.5,alpha=.85)+geom_line(alpha=.65)+
  facet_wrap(~paper,scales='free',ncol=3)+
  scale_color_manual(values=c('#176b9b','#cf5d20'))+
  labs(title='Default tuning, different optimization',
    subtitle='Same outcome samples and penalty; panel axis scales vary',
    x='Observed periods before adoption',y='Time weight',color=NULL,shape=NULL,
    caption='One preselected example per paper. Numerical acceptance and near-flat directions are documented in the report.')+
  theme_minimal(base_size=11)+theme(legend.position='top',panel.grid.minor=element_blank())
print(p)
ggsave('default-vs-optimized.pdf',p,width=12,height=13)
ggsave('default-vs-optimized.png',p,width=12,height=13,dpi=150)
rep_summary <- comparison[comparison$representative_default,c('fit_id','paper_default','T0_default',
  'gap_default','gap_qp','accepted_qp','support_1e6_default','support_1e6_qp',
  'l1_from_default_qp','recent_mass_change','objective_reduction')]
write.csv(rep_summary,'representative-summary.csv',row.names=FALSE)
kable(rep_summary[,-1],digits=5)
```

## Full vectors and reproducibility

The original default results and the main collection remain unchanged. All
comparisons are exploratory. The two zero-noise cases are logged without
substitute weights. Each cohort/outcome is a separate problem; the many fits
within papers are not independent studies.

- [All methods' time weights, including raw QP values](time-weights.csv)
- [Fit diagnostics and objective gaps](fits.csv)
- [Default versus optimized comparisons](default-vs-optimized.csv)
- [Support, runs and span at four cutoffs](support-thresholds.csv)
- [QP initializer sensitivity](initialization-sensitivity.csv)
- [Representative comparisons](representative-summary.csv)
- [Method summary](method-summary.csv)
- [Profile PDF](default-vs-optimized.pdf)
- [Failures](failures.csv)
- [Implementation](implementation.json)
- [Executable Quarto source](analysis.qmd)

Run this document with the prepared matrices from the frozen default analysis and
the pinned author package. Every numerical solver, diagnostic and plot appears
inline. No fitted vectors are imported as the basis of the new analysis; the
previous defaults are read only for an independent equality check.