1  Gsynth

This chapter demonstrates the default gsynth estimator—estimator = "gsynth", (Xu 2017)— which estimates latent factors from the control group only and then imputes counterfactuals for treated units. We cover three settings: block DID, staggered DID, and unbalanced panels, with a dedicated section on inference and one on the placebo test.

For plot customization, see Chapter 3. For IFE-EM and matrix completion (the soft-deprecated paths), see Chapter 2.

ImportantConduct diagnostics

GSC imputes treated counterfactuals by projecting treated-unit loadings onto a factor space estimated from controls. The estimator can return tight-looking confidence intervals even when (a) the treated loadings sit outside the convex hull of control loadings, so the imputation is an extrapolation, or (b) the model overfits, and the parallel-trends-on-factors assumption fails. A point estimate plus a CI is not enough. Two checks should accompany every gsynth fit:

• Loading-overlap (plot(out, type = "loading.overlap")): are treated points inside the shaded control hull? Outside-hull treated units flag extrapolation; report results with that caveat or drop those units.
• Placebo test (placeboTest = TRUE, placebo.period = c(...)): does the model recover near-zero effects on held-out pre-treatment periods? A rejected placebo is evidence against the valid assumptions, and the headline ATT should not be read at face value.

Both diagnostics are demonstrated below for the simulated data and again for the staggered turnout application.

library(gsynth)
library(panelView)
data(gsynth)

1.1 Block DID

We start with a simulated dataset: 5 treated units, 45 control units, and 30 time periods. Treatment begins at Period 21 for all treated units.

head(simdata)

1.1.1 Visualizing data

Before estimation, visualize the data structure with panelView.

panelview(Y ~ D, data = simdata, index = c("id", "time"),
          pre.post = TRUE)

Outcome trajectories by treatment status:

panelview(Y ~ D, data = simdata, index = c("id", "time"),
          type = "outcome")

1.1.2 Estimation

We estimate the model using the outcome \(Y\), the treatment indicator \(D\), and two covariates \(X_1\) and \(X_2\).

The force option specifies fixed effects: "none", "unit", "time", or "two-way". Cross-validation (CV = TRUE) selects the number of factors within the range specified by c(0, 5). From v1.5.0 the CV default is rolling-window with k = 20 folds, cv.prop = 0.1, cv.buffer = 1. To reproduce the pre-v1.5.0 block-CV behavior in existing code, set cv.method = "block" and k = 5 explicitly.

out <- gsynth(Y ~ D + X1 + X2, data = simdata,
              index = c("id", "time"), force = "two-way",
              CV = TRUE, r = c(0, 5), se = TRUE,
              inference = "parametric", nboots = 1000,
              parallel = TRUE, cores = 16)

The output reports sigma2 (error variance), IC (information criterion), and MSPE (mean squared prediction error). Cross-validation selects the r that minimizes MSPE.

out$est.att
#>             ATT      S.E.   CI.lower    CI.upper      p.value count
#> -19  0.05868370 0.5781415 -1.0744527  1.19182014 9.191503e-01     5
#> -18 -0.11741144 0.5208471 -1.1382531  0.90343020 8.216495e-01     5
#> -17 -0.89029173 0.5056829 -1.8814120  0.10082857 7.831069e-02     5
#> -16  1.29461792 0.4913790  0.3315329  2.25770298 8.422087e-03     5
#> -15 -1.20349647 0.5588965 -2.2989134 -0.10807955 3.129163e-02     5
#> -14  2.21991721 0.7193515  0.8100141  3.62982034 2.028704e-03     5
#> -13 -0.18468804 0.7337084 -1.6227300  1.25335396 8.012586e-01     5
#> -12  0.87413414 0.5188940 -0.1428793  1.89114761 9.206377e-02     5
#> -11 -0.04824208 0.6541273 -1.3303081  1.23382393 9.412091e-01     5
#> -10  0.84503247 0.5697891 -0.2717337  1.96179863 1.380580e-01     5
#> -9   0.49754000 0.6708371 -0.8172766  1.81235665 4.582871e-01     5
#> -8  -0.38001372 0.6095868 -1.5747819  0.81475442 5.330246e-01     5
#> -7  -1.03492639 0.5194697 -2.0530682 -0.01678457 4.634090e-02     5
#> -6  -0.81531903 0.7315628 -2.2491557  0.61851764 2.650692e-01     5
#> -5  -1.13328625 0.6272056 -2.3625867  0.09601416 7.078075e-02     5
#> -4  -0.17733004 0.5792112 -1.3125631  0.95790306 7.594845e-01     5
#> -3   0.17670638 0.8458684 -1.4811652  1.83457797 8.345222e-01     5
#> -2  -0.14376360 0.6291828 -1.3769393  1.08941210 8.192633e-01     5
#> -1  -1.19696676 0.5186449 -2.2134921 -0.18044142 2.100618e-02     5
#> 0    1.35910371 0.8934824 -0.3920896  3.11029699 1.282270e-01     5
#> 1    0.34558474 0.6291077 -0.8874436  1.57861309 5.827822e-01     5
#> 2    0.15709192 0.8194300 -1.4489614  1.76314520 8.479703e-01     5
#> 3    3.64127124 0.5989598  2.4673316  4.81521082 1.206894e-09     5
#> 4    2.74975180 0.7548782  1.2702176  4.22928596 2.698532e-04     5
#> 5    4.45569195 0.7889446  2.9093890  6.00199487 1.626449e-08     5
#> 6    5.67231032 0.7678700  4.1673129  7.17730779 1.501022e-13     5
#> 7    8.17038502 0.5937761  7.0066052  9.33416483 0.000000e+00     5
#> 8    6.57084917 0.9056343  4.7958386  8.34585975 4.001244e-13     5
#> 9    9.11624492 0.5342584  8.0691177 10.16337217 0.000000e+00     5
#> 10   9.96204576 0.5057747  8.9707456 10.95334597 0.000000e+00     5
out$est.avg
#>       ATT.avg      S.E. CI.lower CI.upper p.value
#> [1,] 5.084123 0.3336222 4.430235  5.73801       0
out$est.beta
#>        Coef       S.E. CI.lower CI.upper p.value
#> X1 1.454447 0.04319495 1.369787 1.539108       0
#> X2 3.507197 0.04121916 3.426409 3.587985       0

1.1.3 Plot results

The default plot() produces a “gap” plot showing the estimated ATT by period. The true effects in simdata go from 1 to 10 (plus noise) from period 21 to 30.

plot(out)

For counterfactual plots, factor / loading plots, and customization (axes, legends, modern vs legacy visual recipe, highlight API), see Chapter 3.

1.1.4 Factors and loadings

We fit with two factors (r=2) to inspect the estimated factor structure directly:

out_r2 <- gsynth(Y ~ D + X1 + X2, data = simdata,
                 index = c("id", "time"), force = "two-way",
                 r = 2, CV = FALSE, se = FALSE,
                 parallel = FALSE)

type = "factors" plots the estimated time-varying factors (one curve per factor), and type = "loadings" shows the distribution of unit-level loadings on each factor:

plot(out_r2, type = "factors", xlab = "Time")

plot(out_r2, type = "loadings")

1.1.5 Loading-overlap diagnostic

The generalized synthetic control method imputes each treated unit’s counterfactual by projecting its loadings onto a factor space estimated from controls. When the treated-unit loadings sit outside the convex hull of control-unit loadings, the imputed counterfactual is an extrapolation rather than a weighted combination of observed control trajectories.

plot(out, type = "loading.overlap") shows treated and control loadings in the first two factor dimensions, with the convex hull of control loadings shaded. Use it as an at-fit diagnostic, paired with the placebo test below.

plot(out_r2, type = "loading.overlap")

When all treated points fall inside the shaded control hull, the GSC counterfactual is supported by interpolation across observed control trajectories. When treated points fall outside the hull, the counterfactual extrapolates beyond the observed support, and the placebo test below should be checked carefully.

1.2 Inference

gsynth() supports three inference methods via the inference argument; the default depends on the estimator:

inference =	Method	When to use	Default for
"parametric"	Parametric bootstrap with AR-aware residual draws (Xu (2017); AR-preserving fix in fect v2.3.0)	Small \(N_{tr}\); finite-sample setting	estimator = "gsynth"
"bootstrap" (alias "nonparametric")	Block bootstrap at unit level (or cluster level if cl set)	Large \(N\), \(N_{tr}\) growing with \(N\); super-population framework	estimator = "ife" and "mc"
"jackknife"	Leave-one-unit-out jackknife	Cluster-robust under leave-one-unit-out resampling; computationally cheap when \(N_{co}\) is moderate	—

1.2.1 Confidence-interval method (ci.method)

New in v1.5.0. The ci.method argument controls how the confidence intervals stored in the est.* slots (e.g. out$est.att) are formed:

• ci.method = "normal" (default): Wald intervals, \(\hat{\theta} \pm z_{1-\alpha/2}\,\widehat{\mathrm{SE}}\). Symmetric around the point estimate.
• ci.method = "basic": reflected-pivot bootstrap interval (Davison and Hinkley (1997), Sec. 5.2.1). Asymmetric; reads tail quantiles of the bootstrap distribution. Requires inference = "bootstrap" or "jackknife". With \(B < 1000\) replications the tail-quantile endpoints can be unstable; fect emits a warning and recommends nboots = 1000 for publication-grade CIs.

For percentile, BC, and BCa intervals on the alternative estimands (att.cumu, aptt, log.att), call estimand(fit, type, ci.method = ...) post-fit.

1.2.2 Number of bootstrap replicates (nboots)

The default is nboots = 200, sufficient for SEs and Wald ("normal") CIs. For tail-quantile-based intervals (ci.method = "basic", "percentile", "bc", "bca"), bump to nboots = 1000+.

1.2.3 Parallel computing

parallel = TRUE (default) speeds up both the bootstrap and the cross-validation. Set cores explicitly when running on a shared machine; if cores = NULL, gsynth selects min(8, parallel::detectCores() - 2).

out_fast <- gsynth(Y ~ D + X1 + X2, data = simdata,
                   index = c("id", "time"), force = "two-way",
                   CV = TRUE, r = c(0, 5), se = TRUE,
                   inference = "parametric", nboots = 1000,
                   parallel = TRUE, cores = 16)
print(out_fast)
#> Call:
#> gsynth(formula = Y ~ D + X1 + X2, data = simdata, index = c("id", 
#>     "time"), force = "two-way", r = c(0, 5), CV = TRUE, se = TRUE, 
#>     nboots = 1000, inference = "parametric", parallel = TRUE, 
#>     cores = 16, vartype = "parametric")
#> 
#> Average Treatment Effect on the Treated:
#>      ATT.avg   S.E. CI.lower CI.upper p.value
#> [1,]   5.544 0.2609    5.032    6.055       0
#> 
#>    ~ by Period (including Pre-treatment Periods):
#>           ATT   S.E. CI.lower CI.upper   p.value count
#> -19  0.392160 0.5578  -0.7012   1.4855 4.821e-01     5
#> -18  0.276548 0.4431  -0.5920   1.1451 5.326e-01     5
#> -17 -0.275393 0.5190  -1.2926   0.7418 5.957e-01     5
#> -16  0.441201 0.4447  -0.4305   1.3129 3.212e-01     5
#> -15 -0.889595 0.5115  -1.8921   0.1129 8.199e-02     5
#> -14  0.593891 0.4603  -0.3084   1.4961 1.970e-01     5
#> -13  0.528749 0.4090  -0.2729   1.3304 1.961e-01     5
#> -12  0.171569 0.5432  -0.8931   1.2362 7.521e-01     5
#> -11  0.610832 0.4920  -0.3534   1.5751 2.144e-01     5
#> -10  0.170597 0.4553  -0.7219   1.0630 7.079e-01     5
#> -9  -0.271892 0.5805  -1.4097   0.8659 6.395e-01     5
#> -8   0.094843 0.4815  -0.8488   1.0385 8.438e-01     5
#> -7  -0.651976 0.5506  -1.7312   0.4272 2.364e-01     5
#> -6   0.573686 0.4677  -0.3430   1.4904 2.200e-01     5
#> -5  -0.469686 0.5180  -1.4850   0.5456 3.646e-01     5
#> -4  -0.077766 0.5715  -1.1978   1.0423 8.918e-01     5
#> -3  -0.141785 0.5439  -1.2078   0.9242 7.943e-01     5
#> -2  -0.157100 0.4055  -0.9518   0.6376 6.984e-01     5
#> -1  -0.915575 0.5337  -1.9616   0.1304 8.624e-02     5
#> 0   -0.003309 0.3431  -0.6757   0.6691 9.923e-01     5
#> 1    1.235962 0.7414  -0.2171   2.6890 9.548e-02     5
#> 2    1.630264 0.5744   0.5044   2.7561 4.538e-03     5
#> 3    2.712178 0.5787   1.5779   3.8464 2.779e-06     5
#> 4    3.466758 0.7096   2.0759   4.8576 1.033e-06     5
#> 5    5.740132 0.5530   4.6562   6.8241 0.000e+00     5
#> 6    5.280035 0.5486   4.2048   6.3553 0.000e+00     5
#> 7    8.436485 0.4748   7.5059   9.3671 0.000e+00     5
#> 8    7.839902 0.6418   6.5819   9.0979 0.000e+00     5
#> 9    9.455115 0.5411   8.3946  10.5157 0.000e+00     5
#> 10   9.638509 0.4656   8.7259  10.5512 0.000e+00     5
#> 
#> Coefficients for the Covariates:
#>     Coef    S.E. CI.lower CI.upper p.value
#> X1 1.022 0.02911   0.9648    1.079       0
#> X2 3.053 0.03065   2.9929    3.113       0

1.2.4 Comparing CI methods on the same fit

out_basic <- gsynth(Y ~ D + X1 + X2, data = simdata,
                    index = c("id", "time"), force = "two-way",
                    CV = TRUE, r = c(0, 5), se = TRUE,
                    inference = "bootstrap", nboots = 1000,
                    ci.method = "basic",
                    parallel = TRUE, cores = 16)
head(out_basic$est.att)
#>            ATT      S.E.   CI.lower   CI.upper p.value count
#> -19  0.0586837 0.7005020 -1.2622884  1.4069323   0.906     5
#> -18 -0.1174114 0.3997722 -0.8984315  0.6509710   0.738     5
#> -17 -0.8902917 0.2465447 -1.3557802 -0.4014071   0.000     5
#> -16  1.2946179 0.2979121  0.7011996  1.8508511   0.000     5
#> -15 -1.2034965 0.5228902 -2.2618413 -0.2887846   0.022     5
#> -14  2.2199172 0.5887002  1.0820618  3.2791318   0.000     5

The point estimate and SE are identical to a ci.method = "normal" fit on the same bootstrap replicates; only the CI columns differ. Asymmetry of the basic interval relative to the Wald interval signals skewness in the bootstrap distribution.

1.3 Placebo test

A placebo test re-fits the model after pretending some pre-treatment periods are part of the treated window, then asks whether the estimated effect on those held-out periods is close to zero. If the model is well-specified and the parallel-trends assumption holds, the placebo effect should be near zero. If it is not, the headline ATT estimate is suspect.

In v1.5.0, gsynth() accepts two new arguments for this:

• placeboTest = TRUE — turn the placebo test on.
• placebo.period = c(start, end) — which event-time periods to use as the placebo block. For example, c(-2, 0) treats event-times \(-2\) through \(0\) as placebo.

The default plot() highlights the placebo region automatically when placeboTest = TRUE was set at fit time.

out_pb <- gsynth(Y ~ D + X1 + X2, data = simdata,
                 index = c("id", "time"), force = "two-way",
                 r = 2, CV = FALSE, se = TRUE,
                 inference = "bootstrap", nboots = 500,
                 placeboTest = TRUE, placebo.period = c(-2, 0),
                 parallel = TRUE, cores = 16)

plot(out_pb, main = "Placebo test on simdata")

The plot has two regions worth checking:

• Pre-placebo (event time \(< -2\)): the estimated effect should fluctuate around zero. Persistent deviation here means the model is fitting noise as signal, not an overlap problem.
• Placebo region (event times \(-2\) through \(0\)): the confidence interval should cover zero. If it does, the test fails to reject pre-treatment zero, which is consistent with the parallel-trends assumption.

For numeric output by event time, the slot out_pb$est.placebo holds the placebo-region inference (point estimate, SE, CI), and out_pb$est.att holds the ATT estimates outside the placebo block.

out_pb$est.placebo
#>                      Coef      S.E.   CI.lower    CI.upper    p.value
#> Placebo effect -0.4930715 0.2069109 -0.8986094 -0.08753362 0.01717155
#>                CI.lower(90%) CI.upper(90%)
#> Placebo effect    -0.8334096    -0.1527334

A note on language: a placebo test that fails to reject pre-treatment zero is consistent with the assumption; it does not prove the assumption holds. Conversely, a placebo test that rejects pre-treatment zero is evidence against the model’s pre-treatment fit, and the post-treatment ATT should then be reported with that caveat.

The placebo test comes back in Chapter 2 applied to the IFE-EM and matrix-completion estimators.

1.4 Staggered DID

The turnout dataset examines the effect of Election Day Registration (EDR) reforms on voter turnout in the United States, where treatment kicks in at different times.

names(turnout)
#> [1] "abb"            "year"           "turnout"        "policy_edr"    
#> [5] "policy_mail_in" "policy_motor"

1.4.1 Data structure

The balanced panel has 9 treated units with staggered treatment timing.

panelview(turnout ~ policy_edr, data = turnout,
          index = c("abb", "year"), pre.post = TRUE,
          by.timing = TRUE)

Outcome trajectories by treatment-status group:

panelview(turnout ~ policy_edr, data = turnout,
          index = c("abb", "year"), type = "outcome",
          main = "EDR Reform and Turnout",
          by.group = TRUE)

1.4.2 Estimation without factors

When no factor is assumed, the estimates are close to difference-in-differences:

out0 <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
               data = turnout, index = c("abb", "year"),
               se = TRUE, inference = "parametric",
               r = 0, CV = FALSE, force = "two-way",
               nboots = 1000, seed = 02139)

plot(out0, type = "gap", xlim = c(-10, 5), ylim = c(-30, 20))

1.4.3 Estimation with factors

Cross-validation selects the number of factors:

out <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
              data = turnout, index = c("abb", "year"),
              se = TRUE, inference = "parametric",
              r = c(0, 5), CV = TRUE, force = "two-way",
              nboots = 1000, seed = 02139)

1.4.4 Implied weights

out$wgt.implied (N_co × N_tr) stores the implied weights of all control units for each treated unit. Unlike classical synthetic control, the weights can be both positive and negative.

dim(out$wgt.implied)
#> [1] 38  9
sort(out$wgt.implied[, 8])
#>  [1] -0.448332996 -0.440735855 -0.392691119 -0.357515853 -0.351555588
#>  [6] -0.315036478 -0.307735054 -0.262840490 -0.249508765 -0.241817058
#> [11] -0.196655070 -0.145968922 -0.100494444 -0.097927058 -0.068379158
#> [16] -0.061107155 -0.027163955 -0.021590866  0.002553607  0.024038653
#> [21]  0.059699609  0.076882137  0.108178439  0.112349196  0.112475584
#> [26]  0.117607833  0.180530853  0.197468762  0.214643936  0.233002776
#> [31]  0.260495419  0.275884407  0.289004505  0.292635343  0.340999767
#> [36]  0.362079663  0.394710893  0.431814501

1.4.5 Plot results

plot(out, type = "gap", xlim = c(-10, 5), ylim = c(-15, 15))

State-specific gap (Wisconsin):

plot(out, type = "gap", id = "WI", main = "Wisconsin")

1.4.6 Factors and loadings

The estimated time-varying factors and unit-level loadings:

plot(out, type = "factors", xlab = "Year")

plot(out, type = "loadings")

Loading-overlap check on the turnout fit, in the first two factor dimensions:

plot(out, type = "loading.overlap")

Additional plot types and customization (raw, counterfactual, modern theme, highlight API) are in Chapter 3.

1.5 Unbalanced panels

gsynth can accommodate unbalanced panels. We demonstrate by randomly removing observations from the turnout dataset.

set.seed(123456)
turnout.ub <- turnout[-c(which(turnout$abb == "WY")[1:15],
                         sample(1:nrow(turnout), 50,
                                replace = FALSE)), ]

Visualize the data structure. White cells represent missing values.

panelview(turnout ~ policy_edr + policy_mail_in + policy_motor,
          data = turnout.ub, index = c("abb", "year"),
          pre.post = TRUE)

1.5.1 Estimation

out <- gsynth(turnout ~ policy_edr + policy_mail_in + policy_motor,
              data = turnout.ub, index = c("abb", "year"),
              se = TRUE, inference = "parametric",
              r = c(0, 5), CV = TRUE, force = "two-way",
              parallel = TRUE, min.T0 = 8,
              nboots = 1000, seed = 02139)

On unbalanced panels, the parametric bootstrap uses a unit-level Gaussian approximation for residual perturbations (introduced in fect v2.3.0); residuals are drawn from a Gaussian fitted to each unit’s residual structure rather than from the empirical pool, naturally accommodating missing cells. This is automatic — no API change for the user. Balanced panels see no change in behavior.

1.5.2 Missing-data plot

The "missing" plot shows the data used in estimation. Gray cells with red dots indicate observations removed due to insufficient pre-treatment periods. The matrix is available in out$obs.missing (0 = missing, 1 = control, 2 = treated pre, 3 = treated post, 4 = treated removed).

plot(out, type = "missing")

Subset to treated units only:

plot(out, type = "missing", xlab = "Year", ylab = "State",
     main = "Treatment Status", id = out$id.tr,
     xlim = c(1920, 2012), axis.adjust = TRUE)

Gap plot with the unbalanced panel:

plot(out, type = "gap", xlim = c(-10, 5), ylim = c(-15, 15))

1.5.3 Notes on performance

1. Unbalanced panels take significantly more time than balanced panels (often 100:1 or higher) due to the EM algorithm for missing values.
2. Adding covariates also slows the algorithm because the IFE/MC model takes more iterations to converge.
3. Setting min.T0 to a positive number helps — the algorithm discards treated units with too few pre-treatment periods. A larger T_0 reduces bias and extrapolation risk. For cross-validation, min.T0 must be at least r.max + 2.

Athey, Susan, Mohsen Bayati, Nikolay Doudchenko, Guido Imbens, and Khashayar Khosravi. 2021. “Matrix Completion Methods for Causal Panel Data Models.” Journal of the American Statistical Association forthcoming.

Chen, Jiafeng, and Jonathan Roth. 2024. “Logs with Zeros? Some Problems and Solutions.” Quarterly Journal of Economics.

Davison, A. C., and D. V. Hinkley. 1997. Bootstrap Methods and Their Application. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press.

Gobillon, Laurent, and Thierry Magnac. 2016. “Regional Policy Evaluation: Interactive Fixed Effects and Synthetic Controls.” Review of Economics and Statistics 98 (3): 535–51.

Xu, Yiqing. 2017. “Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models.” Political Analysis 25 (1): 57–76.