Causal Inference VII(b) — Synthetic Difference-in-Differences (R companion)¶

From-scratch SDID via quadprog, with fixest (DiD)¶

The Arkhangelsky et al. synthdid package is GitHub-only, so this companion implements synthetic DiD from scratch in R, using quadprog to solve the constrained (simplex + ridge) weight problems exactly — matching the Python notebook to the dollar. It compares SDID to DiD (fixest) and a from-scratch synthetic control on the Proposition 99 data.

1. The three estimators — DiD, SC, and SDID¶

The unit and time weights are constrained least-squares problems (non-negative, summing to one, with ridge regularization for SDID's unit weights), solved with quadprog::solve.QP. DiD (equal weights, fixest) overstates the effect; synthetic control and synthetic DiD are more credible, SDID reproducing the published ≈ −15.6.

In [1]:
suppressMessages({library(fixest); library(quadprog)})
d<-read.csv("smoking.csv")
piv<-reshape(d[,c("state","year","cigsale")], idvar="year", timevar="state", direction="wide")
yrs<-piv$year; piv$year<-NULL; colnames(piv)<-sub("cigsale.","",colnames(piv)); M<-as.matrix(piv)
tr<-"California"; Y1<-M[,tr]; Y0<-M[,colnames(M)!=tr]; donors<-colnames(Y0)
pre<-yrs<=1988; post<-yrs>=1989; Npre<-sum(pre); Npost<-sum(post)
# QP: min ||A w - b||^2 s.t. sum(w)=1, w>=0 (+ ridge)
qp_simplex<-function(A,b,ridge=0){ n<-ncol(A); D<-t(A)%*%A+diag(ridge,n)+diag(1e-8,n)
  solve.QP(D, t(A)%*%b, cbind(rep(1,n),diag(n)), c(1,rep(0,n)), meq=1)$solution }
# DiD via fixest
dl<-do.call(rbind, lapply(colnames(M), function(s) data.frame(state=s, year=yrs, y=M[,s], ca=as.integer(s==tr), post=as.integer(yrs>=1989))))
did<-coef(feols(y~ca:post | state+year, data=dl))["ca:post"]
# synthetic control (unit weights, exact pre-fit)
wsc<-qp_simplex(Y0[pre,], Y1[pre]); sc<-mean(Y1[post]-Y0[post,]%*%wsc)
# synthetic DiD: unit weights (centered pre -> absorbs intercept; ridge) + time weights
zeta<-(Npost)^0.25*sd(apply(Y0[pre,],2,diff))
sdid_est<-function(Y0m){ n<-ncol(Y0m)
  Ap<-sweep(Y0m[pre,],2,colMeans(Y0m[pre,])); bp<-Y1[pre]-mean(Y1[pre])        # center -> level-shift intercept
  om<-qp_simplex(Ap, bp, ridge=zeta^2*Npre)
  At<-t(Y0m[pre,]); lam<-qp_simplex(sweep(At,2,colMeans(At)), colMeans(Y0m[post,])-mean(colMeans(Y0m[post,])))
  (mean(Y1[post])-sum(lam*Y1[pre]))-(sum(colMeans(Y0m[post,])*om)-sum(lam*(Y0m[pre,]%*%om))) }
om<-{Ap<-sweep(Y0[pre,],2,colMeans(Y0[pre,])); qp_simplex(Ap,Y1[pre]-mean(Y1[pre]),ridge=zeta^2*Npre)}
lam<-{At<-t(Y0[pre,]); qp_simplex(sweep(At,2,colMeans(At)), colMeans(Y0[post,])-mean(colMeans(Y0[post,])))}
sdid<-sdid_est(Y0)
cat(sprintf("DiD (fixest)          = %.2f\n", did))
cat(sprintf("Synthetic Control     = %.2f\n", sc))
cat(sprintf("Synthetic DiD (SDID)  = %.2f   (Python -15.60; Arkhangelsky ~ -15.6)\n", sdid))
options(repr.plot.width=7.5, repr.plot.height=4)
bp<-barplot(c(DiD=did,`Synth Control`=sc,`Synthetic DiD`=sdid), col=c("#a0aec0","#2b6cb0","#2f855a"), ylab="estimated effect (packs)", main="Prop 99: DiD vs SC vs Synthetic DiD"); abline(h=0)
text(bp, c(did,sc,sdid)-1.5, sprintf("%.1f",c(did,sc,sdid)), col="white", font=2)

Warning message:
"package 'fixest' was built under R version 4.6.1"

DiD (fixest)          = -27.35

Synthetic Control     = -19.51

Synthetic DiD (SDID)  = -15.61   (Python -15.60; Arkhangelsky ~ -15.6)
No description has been provided for this image

2. Weights and jackknife inference¶

The unit weights (a sparse set of donor states) and time weights (emphasizing the late-1980s), plus a leave-one-control-out jackknife standard error. With the exact QP solver the jackknife is stable and matches the Python notebook (SE ≈ 2.4; the effect is significant).

In [2]:
options(repr.plot.width=13, repr.plot.height=4.2); par(mfrow=c(1,2))
top<-sort(setNames(om,donors), decreasing=TRUE)[1:8]
barplot(rev(top), horiz=TRUE, las=1, col="#2f855a", xlab="unit weight", main="SDID unit weights (top states)")
barplot(lam, names.arg=yrs[pre], col="#6b46c1", xlab="pre-year", ylab="time weight", main="SDID time weights")
par(mfrow=c(1,1))
jk<-sapply(1:ncol(Y0), function(j) sdid_est(Y0[,-j,drop=FALSE])); Nc<-ncol(Y0)
se<-sqrt((Nc-1)/Nc*sum((jk-mean(jk))^2))
cat(sprintf("SDID = %.2f;  jackknife SE = %.2f;  95%% CI = [%.2f, %.2f]  (significant)\n", sdid, se, sdid-1.96*se, sdid+1.96*se))

SDID = -15.61;  jackknife SE = 2.42;  95% CI = [-20.34, -10.87]  (significant)
No description has been provided for this image

3. Summary¶

Using quadprog to solve the weight problems exactly, R reproduced the Python notebook to the dollar: DiD (−27.4), synthetic control (−19.5), and synthetic DiD (−15.6, matching Arkhangelsky et al.), with a stable jackknife CI. SDID's twin reweightings — unit weights for trend, time weights for the relevant pre-period, plus the level-shift double difference — make it more robust than either parent.

For production use, the authoritative implementation is the synthdid package (GitHub); Synth/tidysynth cover classic synthetic control and fixest covers DiD. Recommendation: use synthetic DiD as a default for panel comparative-case studies, reported beside SC and DiD. Cross-links: SDID unifies Synthetic Control (subsection 7) and Difference-in-Differences (subsection 6); its doubly-robust flavour echoes AIPW/DML/TMLE. This completes the depth of the Synthetic Control subsection.