6 Gsynth Program

This chapter demonstrates the generalized synthetic control method, or Gsynth, proposed in Xu (2017) [Paper]. R script used in this chapter can be downloaded here.

The gsynth package as a standalone package still exists, but algorithmically is a wrapper of the fect package. Click here if you want to use the old syntax.

Recommended changes in v2.3.0 for Gsynth users

Two new defaults are strongly recommended for the Synth setting:

1. Use rolling-window CV for rank selection. This is now the default. Gsynth applications often use macro panels, where residual shocks are persistent over time. Before v2.3.0, block CV could leak information across the train/test boundary and over-select r, often choosing r.max. Rolling CV (cv.method = "rolling") prevents this future-side leakage by design.

2. Bound treated loadings with loading.bound = "simplex". This restricts each treated unit’s factor loading to a non-negative weighted average of the control loadings. In other words, the imputed counterfactual is kept inside the support of the donor pool, matching the usual synthetic-control intuition. The default, loading.bound = "none", preserves the original Xu (2017) behavior: an unconstrained least-squares projection. The simplex bound is opt-in, but it is generally preferable, especially when the donor pool is small or overlap with the donor pool is uncertain. See the Bounded factor loadings section below and the companion loading.overlap diagnostic.

Gsynth was originally implemented in the gsynth package but has now been fully integrated into the fect package. Gsynth (method = "gsynth") and FEct/IFEct/MC (method = "fe"/"ife"/"mc") are different in the following ways:

Gsynth is designed to handle block and staggered DID settings without treatment reversal, whereas other methods allow for treatment reversal under the assumption of limited carryover effects.
Gsynth is particularly suited for the Synth setting, where the number of treated units is small, including scenarios with only one treated unit. By setting vartype = "parametric", we can use a two-stage parametric bootstrapping procedure to produce uncertainty estimates. In contrast, other methods rely on large samples, particularly a large number of treated units, to obtain reliable standard errors and confidence intervals using "bootstrap" or "jackknife".
Compared with IFEct (method = "ife"), Gsynth does not rely on pre-treatment data from the treated units to estimate time components (e.g., factors). Hence, time.component.from = "nevertreated". This approach speeds up computation, improves stability, and is more suitable for predictive inference based on row exchangeability.

Therefore, we recommend setting method = "gsynth" in fect for the synthetic control setting, where the treatment does not reverse (or is coded accordingly) and the number of treated units is small.

Key Equivalence

fect(..., method = "gsynth") is equivalent to fect(..., method = "ife", time.component.from = "nevertreated").

Gsynth is a procedure based on an IFE model estimated using only never-treated controls. The gsynth package was developed in the synthetic control setting, where a small number of treated units are compared with never-treated controls to learn latent time effects.

Use method = "gsynth" (or equivalently time.component.from = "nevertreated") when:

The treatment does not reverse (block or staggered adoption only).
The number of treated units is small, including $N_{tr} = 1$.
You want parametric bootstrap (vartype = "parametric") for uncertainty estimation — this is particularly useful with few treated units.
You do not want pre-treatment data from treated units to influence time component (i.e., factor) estimation.

In contrast, the default regime (time.component.from = "notyettreated", used in Chapter 2) accommodates treatment reversal and uses EM to estimate time components (time fixed effects, latent factors, and temporal dynamics) from all not-yet-treated observations.

We will use two datasets, sim_gsynth and turnout, to perform analyses in block and staggered DID settings. First, we load the two datasets: sim_gsynth and turnout:

data(sim_gsynth)
data(turnout)
ls()
#> [1] "sim_gsynth" "simdata"    "turnout"

6.1 Basic Usage

6.1.1 Simulated Data

We start with the first example, sim_gsynth, a simulated dataset described in Xu (2017).

There are 5 treated units, 45 control units, and 30 time periods. The treatment kicks at Period 21 for all treated units, hence a block DID setup.

head(sim_gsynth)
#>    id time         Y D        X1         X2 eff      error mu       alpha
#> 1 101    1  6.210998 0 0.3776736 -0.1732470   0  0.2982276  5 -0.06052364
#> 2 101    2  4.027106 0 1.7332009 -0.4945009   0  0.6365697  5 -0.06052364
#> 3 101    3  8.877187 0 1.8580159  0.4984432   0 -0.4837806  5 -0.06052364
#> 4 101    4 11.515346 0 1.3943369  1.1272713   0  0.5168620  5 -0.06052364
#> 5 101    5  5.971526 0 2.3636963 -0.1535215   0  0.3689645  5 -0.06052364
#> 6 101    6  8.237905 0 0.5370867  0.8774397   0 -0.2153805  5 -0.06052364
#>           xi          F1          L1           F2         L2
#> 1  1.1313372  0.25331851 -0.04303273  0.005764186 -0.8804667
#> 2 -1.4606401 -0.02854676 -0.04303273  0.385280401 -0.8804667
#> 3  0.7399475 -0.04287046 -0.04303273 -0.370660032 -0.8804667
#> 4  1.9091036  1.36860228 -0.04303273  0.644376549 -0.8804667
#> 5 -1.4438932 -0.22577099 -0.04303273 -0.220486562 -0.8804667
#> 6  0.7017843  1.51647060 -0.04303273  0.331781964 -0.8804667

Before we conduct any statistical analysis, it is helpful to visualize the data structure and spot missing values (if there are any). We can easily do so with the help of the panelView package. The figure below shows that: (1) there are 5 treated units and 45 control units; (2) the treated units start to be treated in period 21; and (3) there are no missing values, which is a rare case.

library(panelView)
#> ## See https://yiqingxu.org/packages/panelview/ for more info.
#> ## Report bugs -> yiqingxu@stanford.edu.
#> 
#> Attaching package: 'panelView'
#> The following object is masked _by_ '.GlobalEnv':
#> 
#>     simdata
#> The following object is masked from 'package:fect':
#> 
#>     simdata
panelview(Y ~ D, data = sim_gsynth,  index = c("id","time"), pre.post = TRUE)

The code chunk below visualizes the trends of outcome variable over the full panel across groups; different colors correspond to unique treatment statuses.

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

6.1.2 Estimation

We estimate the model using only the outcome variable ($Y$), treatment indicator ($D$), covariates ($X_1$) and ($X_2$), and group indicators ($id$) and ($time$).

To implement the Gsynth algorithm proposed in Xu (2017), set method = "gsynth".

system.time(
out <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id","time"), 
            method = "gsynth", force = "two-way", CV = TRUE, r = c(0, 5), 
            se = TRUE, nboots = 1000, vartype = 'parametric', 
            parallel = TRUE, cores = 16))
#> 
#>  +----------------------------------------------------------+
#>  | Parallel computing: using 16 of 14 available cores.      |
#>  |                                                          |
#>  | To change: set cores = <n> in fect().                    |
#>  | Default: min(available - 2, 8).                          |
#>  +----------------------------------------------------------+
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Interactive fixed effects model...
#> Cross-validating ...
#> r = 0; sigma2 = 1.84865; IC = 1.02023; PC = 1.74458; MSPE = 1.87524
#> r = 1; sigma2 = 1.51541; IC = 1.20588; PC = 1.99818; MSPE = 1.68924
#> r = 2; sigma2 = 0.99737; IC = 1.16130; PC = 1.69046; MSPE = 1.30085
#> *
#> r = 3; sigma2 = 0.94664; IC = 1.47216; PC = 1.96215; MSPE = 2.12741
#> r = 4; sigma2 = 0.89411; IC = 1.76745; PC = 2.19241; MSPE = 2.31354
#> r = 5; sigma2 = 0.85060; IC = 2.05928; PC = 2.40964; MSPE = 2.71562
#> 
#>  r* = 2
#> 
#> Parametric Bootstrap (para.error = "empirical")
#> 
Simulating errors ...
#> Can't calculate the F statistic because of insufficient treated units.
#>    user  system elapsed 
#>  18.143   1.003  32.549

The first variable on the right-hand side is the binary treatment indicator, with remaining variables as controls. The index option designates unit and time indicators for fixed effects analysis.

The force option ("none", "unit", "time", and "two-way") specifies the additive component(s) of the fixed effects included in the model. The default option is "two-way" (including additive unit fixed effects). A cross-validation procedure is provided (when CV = TRUE) to select the number of unobserved factors within the interval of r=c(0,5). When cross-validation is switched off, the first element in r will be set as the number of factors. Alternatively, for example, when CV = FALSE, r can take in a numeric value ($r = 0$).

Setting se = TRUE, the algorithm can produce uncertainty measurements. When the number of treated units is small, a parametric bootstrap procedure is preferred (vartype = "parametric"). The default number of parametric bootstrap runs is set to 200. The number of bootstrap runs can be set by nboots.

The algorithm prints out the results automatically. sigma2 stands for the estimated variance of the error term; IC represents the Bayesian Information Criterion; and MPSE is the Mean Squared Prediction Error. The cross-validation procedure selects an r* that minimizes the MSPE.

Users can also use the print function to directly retrieve specific results. Here, we show three examples for the commonly requested estimates.

est.att reports the average treatment effect on the treated (ATT) by period
est.avg shows the ATT averaged over all periods
beta represents the coefficients of the time-varying covariates

print(out)
out$est.att
out$est.avg
out$beta

Treatment effect estimates from each bootstrap run are stored in att.boot, an array whose dimension = (#time periods * #treated * #bootstrap runs).

Parallel computing will speed up both cross-validation and uncertainty estimation significantly. We recommend that users manually set the number of cores using the cores option. If this is not supplied or is NULL, we will automatically select the smaller of 8 and the number of usable system cores minus 2, to prevent excessive use of system resources.

system.time(
out <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id","time"), method = "gsynth", force = "two-way", CV = TRUE, r = c(0, 5), se = TRUE, nboots = 1000, vartype = 'parametric', parallel = TRUE, cores = 16)
)
#> 
#>  +----------------------------------------------------------+
#>  | Parallel computing: using 16 of 14 available cores.      |
#>  |                                                          |
#>  | To change: set cores = <n> in fect().                    |
#>  | Default: min(available - 2, 8).                          |
#>  +----------------------------------------------------------+
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Interactive fixed effects model...
#> Cross-validating ...
#> r = 0; sigma2 = 1.84865; IC = 1.02023; PC = 1.74458; MSPE = 2.28134
#> *
#> r = 1; sigma2 = 1.51541; IC = 1.20588; PC = 1.99818; MSPE = 2.42098
#> r = 2; sigma2 = 0.99737; IC = 1.16130; PC = 1.69046; MSPE = 1.91650
#> *
#> r = 3; sigma2 = 0.94664; IC = 1.47216; PC = 1.96215; MSPE = 2.82113
#> r = 4; sigma2 = 0.89411; IC = 1.76745; PC = 2.19241; MSPE = 3.40386
#> r = 5; sigma2 = 0.85060; IC = 2.05928; PC = 2.40964; MSPE = 4.16308
#> 
#>  r* = 2
#> 
#> Parametric Bootstrap (para.error = "empirical")
#> 
Simulating errors ...
#> Can't calculate the F statistic because of insufficient treated units.
#>    user  system elapsed 
#>  19.124   0.961  33.472

Gsynth in fect also incorporates jackknife method for uncertainty estimates.

out2 <- fect(Y ~ D + X1 + X2, data = sim_gsynth,  index = c("id","time"),
               method = "gsynth", force = "two-way",
               CV = TRUE, r = c(0, 5), se = TRUE,
               vartype = "jackknife",
               parallel = TRUE, cores = 16)

6.1.3 Cross-validation

cv.method = "rolling" is the new default across all CV entry points (DID/TWFE, Gsynth, and fect_mspe). cv.method = "block" is also available; "loo" is exclusive to Gsynth (time.component.from = "nevertreated") because it requires factor estimation to be independent of the held-out treated observations.

	`fect_cv` (DID/TWFE)	CV in Gsynth	`fect_mspe` (post-hoc)
`"rolling"`	Default	Default	Default
`"block"`	Available	Available	Available
`"loo"`	—	Available	—

The legacy values "all_units" (now "block") and "treated_units" (block masking restricted to treated pre-treatment cells) are still accepted but emit a deprecation message; both will be replaced by the unified (cv.method, cv.units) API in v2.4.0. See the cheatsheet for the deprecation plan.

6.1.4 Rolling-window CV

Set cv.method = "rolling" in the main fect() CV dispatcher to use rolling-window CV (the recommended default for serially correlated panels — see the recommended-changes callout at the top of this chapter for the empirical motivation, and Chapter 4 for the design and the side-by-side figure).

fit <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id", "time"),
            method = "gsynth", force = "two-way",
            CV = TRUE, r = c(0, 5),
            cv.method = "rolling",
            cv.buffer = 1, cv.nobs = 3, k = 20, cv.prop = 0.1,
            cv.rule = "1se",
            se = TRUE, vartype = "parametric")
fit$r.cv     # selected r --- use this for downstream inference

CV behavior is identical across method = "ife" (IFE-EM, internal time.component.from = "notyettreated") and method = "gsynth" (GSC, internal time.component.from = "nevertreated"); both populate Y.ct.full at masked positions, so the rolling design scores MSPE uniformly. The full parameter reference for all CV methods lives in the Cross-Validation section of the cheatsheet.

6.1.5 Visualizing Results

By default, the print function produces a gap plot, equivalent to using plot(out, type = "gap"), and visualizes the estimated Average Treatment Effect on the Treated (ATT) by period. For reference, the true effects in gsynth range from 1 to 10 (with some added white noise) and take effect during periods 21 to 30.

Note that the plot() generates a ggplot2 object, which users can save by simply calling print with the object’s name.

a <- plot(out) # by default, type = "gap"
print(a)

By switching on connected, the confidence interval of the ATT estimates will be represented by a shaded area.

plot(out, connected = TRUE)

Moreover, by switching off show.points, ATT estimates can be displayed as a line plot.

plot(out, connected = TRUE, show.points = FALSE)

Aesthetic options for ggplot2 objects are compatible for all fect plots.

For demonstration, we use main, xlim, and ylim to set the plot title and axis labels.

plot(out, type = "gap", ylim = c(-6,12), xlab = "Period", 
     main = "Estimated ATT (Gsynth)")

fect objects can generate eight types of plots for diverse demonstration purposes. The type option includes the following:

"gap": Reports ATT by period (default).
"counterfactual": Generates the observed treated outcome alongside the imputed counterfactual averages. The style can be modified with raw = "all" for linear or raw = "band" for graphic representation of all estimated counterfactual results.
"factors" and "loadings": Plot the estimated factors and loadings, respectively.
"status": Plots treatment status of all observations, similar to panelview(). The time labels are displayed only when axis.lab = "time" is set.
"box": Displays the estimated individual treatment effects with box plot.
"calendar": Shows how the treatment effect evolves over time.
"equiv": Exhibits the average pre-treatment residuals with equivalence confidence intervals.

To visualize the estimated counterfactual outcomes, we use type = "counterfactual". If no argument is assigned to raw, the plot will display the average of the treated and of the estimated counterfactual outcomes. When method = "gsynth" and vartype = "parametric", the 95% confidence intervals around the counterfactual outcomes will be the same as for the event-study ATT (though flipped in shape).

plot(out, type = "counterfactual")

fect offers two options for displaying all estimated counterfactual estimates: raw = "all" for linear and raw = "band" for a graphical (shaded band) representation.

plot(out, type = "counterfactual", raw = "all")

Note that in the plot below, the shaded areas represent the 5-95% quantiles of the treated and estimated counterfactual trajectories, not their uncertainty estimates.

plot(out, type = "counterfactual", raw = "band")

The figure below illustrates the treatment status of all observations. Users can specify the exact values shown on the x- and y-axes using the xticklabels and yticklabels parameters (i.e., xticklabels=c("1","5","10", "15","20", "25", "30")). Setting a nonexistent value to the option removes all numeric labels for that axis (i.e., yticklabels="0"). The status plot has now been incorporated into a standalone package, panelView.

plot(out, type = "status", yticklabels="0", 
     xticklabels=c("5", "10", "15","20", "25", "30") )

The next two figures plot the estimated factors and factor loadings, respectively.

plot(out, type = "loadings")

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

The box plot visualizes the estimated individual treatment effects for the treated units. While these effects are not identified at the individual level, their dispersion provides insight into heterogeneous treatment effects across different time periods and informs model performance. By default, the number of total treated units is labeled on the graph. Note, if the number of treated units is small, the box plot will reduce to a scatter plot, as shown below.

We also provide the xangle and yangle options to allow users to tilt the labels for better display. We will demonstrate the functionality in later examples.

plot(out, type = "box", xlab = "time",
     xticklabels=c("-19", "-15", "-10", "-5","0","5","10") )

If we want to focus on specific periods within the full panel, the xlim option is useful. Here, we narrow the time frame to fifteen periods before and ten periods after the treatment.

plot(out, type = "box", xlim = c(-15, 10), 
     xticklabels=c( "-15", "-10", "-5","0","5","10"))

To explore how the treatment effect evolves over time, we can set type = "calendar".

In the plot below, the points represent the ATTs by calendar time. The blue curve shows a lowess fit of these estimates, with the shaded band indicating the 95% confidence interval. The red horizontal dashed line marks the overall average ATT across all time periods.

plot(out,type = "calendar")

6.2 Staggered DID

The second example examines the impact of Election-Day Registration (EDR) reforms on voter turnout in the United States, with the treatment taking effect at varying times across states. Further details about the turnout dataset can be found in Xu (2017) Section 5.

6.2.1 Data structure

First, we use panelView to visualize the data structure. The following figure shows that (1) we have a balanced panel with 9 treated units and (2) the treatment starts at different time periods.

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

panelView can visualize the outcome variable using colored lines to represent changes in treatment status.

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

6.2.2 Estimation w/o factors

When no factor is assumed, the ATT estimates are close to what we obtain from difference-in-differences.

out0 <- fect(turnout ~ policy_edr + policy_mail_in + policy_motor, 
               data = turnout, index = c("abb","year"), 
               se = TRUE, method = "gsynth",
               r = 0, CV = FALSE, force = "two-way", 
               nboots = 1000, seed = 02139)
#> 
#>  +----------------------------------------------------------+
#>  | Parallel computing: using 8 of 14 available cores.       |
#>  |                                                          |
#>  | To change: set cores = <n> in fect().                    |
#>  | Default: min(available - 2, 8).                          |
#>  +----------------------------------------------------------+
#> Can't calculate the F statistic because of insufficient treated units.

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

6.2.3 Estimation w/ factors

We now allow the algorithm to determine the optimal number of factors that best predict the pre-treatment data with cross-validation. The keep.sims option must be set to TRUE if the user wishes to visualize subsets of the treated units in the “counterfactual” plot.

out_turnout <- fect(turnout ~ policy_edr + policy_mail_in + policy_motor, 
              data = turnout,  index = c("abb","year"), 
              se = TRUE, method = "gsynth", vartype = "parametric",
              r = c(0, 5), CV = TRUE, force = "two-way", 
              nboots = 1000, seed = 02139, keep.sims = TRUE)
#> 
#>  +----------------------------------------------------------+
#>  | Parallel computing: using 8 of 14 available cores.       |
#>  |                                                          |
#>  | To change: set cores = <n> in fect().                    |
#>  | Default: min(available - 2, 8).                          |
#>  +----------------------------------------------------------+
#> Cross-validating ...
#> Criterion: Mean Squared Prediction Error
#> Interactive fixed effects model...
#> Cross-validating ...
#> r = 0; sigma2 = 77.99366; IC = 4.82744; PC = 72.60594; MSPE = 136.08663
#> r = 1; sigma2 = 13.65239; IC = 3.52566; PC = 21.67581; MSPE = 55.91536
#> *
#> r = 2; sigma2 = 8.56312; IC = 3.48518; PC = 19.23841; MSPE = 56.42439
#> r = 3; sigma2 = 6.40268; IC = 3.60547; PC = 18.61783; MSPE = 79.58436
#> r = 4; sigma2 = 4.74273; IC = 3.70145; PC = 16.93707; MSPE = 45.78004
#> *
#> r = 5; sigma2 = 4.02488; IC = 3.91847; PC = 17.05226; MSPE = 69.60277
#> 
#>  r* = 4
#> 
#> Parametric Bootstrap (para.error = "empirical")
#> 
Simulating errors ...
#> Can't calculate the F statistic because of insufficient treated units.

6.2.4 Implied weights

out$wgt.implied ($N_{co}\times N_{tr}$) stores the implied weights of all control units for each treated unit. Different from the synthetic control method, the weights can take both positive and negative values. Below we show the control unit weights for Wisconsin.

dim(out_turnout$wgt.implied)
#> [1] 38  9
sort(out_turnout$wgt.implied[,8])
#>  [1] -0.448353155 -0.440752705 -0.392714893 -0.357504540 -0.351555127
#>  [6] -0.315049295 -0.307745990 -0.262835262 -0.249519316 -0.241829355
#> [11] -0.196666127 -0.145980298 -0.100508796 -0.097925382 -0.068362340
#> [16] -0.061102735 -0.027155244 -0.021581331  0.002552594  0.024047052
#> [21]  0.059700576  0.076876409  0.108179504  0.112348625  0.112475973
#> [26]  0.117609515  0.180545028  0.197481536  0.214657368  0.232999379
#> [31]  0.260495320  0.275882520  0.288995460  0.292648965  0.341018271
#> [36]  0.362079499  0.394714111  0.431834184

6.2.5 Visualizing results

Like we have demonstrated with the multiperiod DID analysis, we also use the gap (default) plot to visualize the ATT.

plot(out_turnout, xlim = c(-10, 5), ylim=c(-10, 10))

For staggered DID, a status plot can be generated after estimation. Here, we use xlab, ylab, and main to modify the axis labels and the graph title, respectively. Additionally, xangle and yangle can tilt the numeric labels to improve readability. Note that the values provided for these options specify the degree of tilt.

plot(out_turnout, type = "status",xlab = "Year", ylab = "State", main = "Treatment Status", 
     xticklabels=c(1920, 1928, 1936, 1944, 1952, 1960, 
                   1968, 1976, 1984, 1992, 2000, 2008), xangle=10)

We can also visualize this using the “counterfactual” plot.

plot(out_turnout, type = "counterfactual")

If we want to zoom in to examine the treatment effect on a specific state, such as Wisconsin, we can do so by setting the id option accordingly. The id option can also accept a vector if one wishes to plot multiple treated units. This plot will only display confidence intervals for the counterfactual estimates if keep.sims = TRUE is set in the estimation step and vartype = "parametric".

plot(out_turnout, type = "counterfactual", id = "WI", main = "Wisconsin")

As we have seen in the last example, due to the small number of treated units, the box plot reduces to a scatter plot.

plot(out_turnout, type = "box", 
     xticklabels=c("-20", "-15", "-10", "-5","0","5","10"))

Here, we visualize the estimated average treatment effects by calendar time.

plot(out_turnout, type = "calendar", ylim = c(-15,15))

Estimated factors and factor loadings are plotted below.

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

plot(out_turnout, type = "loadings")

6.2.6 Cumulative Effects

For staggered adoption without treatment reversals, the cumulative ATT summarizes how the total treatment dose accumulates over relative time — useful for outcomes where each additional treated period adds to a running total (e.g., cumulative voter mobilization). Because out_turnout was estimated with keep.sims = TRUE and vartype = "parametric", the per-cell replicate surface needed for estimand() is available.

cumu.turnout <- estimand(out_turnout, "att.cumu", "event.time")
esplot(cumu.turnout, Period = "event.time",
       Estimate = "estimate", CI.lower = "ci.lo", CI.upper = "ci.hi",
       Count = "n_cells",
       main = "Cumulative Effect of EDR Reform on Turnout",
       ylab = "Cumulative ATT")

See Chapter 3 (Alternative Estimands) for the full estimand() API including parametric variance, log-scale ATT, and APTT.

6.3 Unbalanced panels

Gsynth in fect can accommodate unbalanced panels. To illustrate how it works, we randomly remove 50 observations as well as the first 15 observations of Wyoming from the turnout dataset and then re-estimate the model.

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

Again, before running any regressions, we first plot the data structure and visualize missing values. In the following plot, white cells represent missing values.

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

6.3.1 Estimation

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

6.3.2 Parametric bootstrap on an unbalanced panel

When the number of treated units is small, the parametric bootstrap (vartype = "parametric") is the recommended uncertainty procedure for gsynth. It is available on unbalanced panels too, with a small technical adjustment.

On a balanced panel, the parametric bootstrap resamples whole-$T$ residual vectors from the controls and adds them to the estimated counterfactual. On an unbalanced panel, individual residual vectors may have missing cells, so direct vector resampling is not always feasible. fect therefore replaces the empirical resample with a Gaussian approximation: it estimates the pairwise-complete empirical covariance matrix of the residuals, $\hat\Sigma \in \mathbb{R}^{T \times T}$, and draws each bootstrap perturbation from a multivariate Gaussian $\mathcal{N}(0, \hat\Sigma)$. The off-diagonal entries of $\hat\Sigma$ encode within-unit serial correlation, so the Gaussian draws inherit the same temporal-dependence structure as the data.

The two procedures also condition on different objects. The parametric bootstrap is conditional on the realized factors $\hat F$, loadings $\hat\Lambda$, covariates $X$, and treatment indicators $D$ — it perturbs only the idiosyncratic error $\varepsilon$. The nonparametric vartype = "bootstrap" is unconditional: by resampling units with replacement, each draw redraws $(F_i, \Lambda_i, X_i, D_i)$ together, so its variance also reflects between-unit heterogeneity in those structural components. The two SEs answer different inferential questions and need not coincide even when both are correctly calibrated.

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

Compare the average CI width of the period-by-period ATT under the two SE procedures on the same unbalanced panel:

ci_width <- function(out) mean(out$est.att[, "CI.upper"] -
                                out$est.att[, "CI.lower"], na.rm = TRUE)
data.frame(
  vartype  = c("bootstrap", "parametric"),
  CI_width = c(ci_width(out_ub), ci_width(out_ub_param))
)
#>      vartype CI_width
#> 1  bootstrap 6.043730
#> 2 parametric 9.713193

Both should be of the same order of magnitude when both correctly capture the panel’s clustering and serial-correlation structure.

6.3.3 Visualizing results

The status plot visualizes the regression data and differs from the panelView plot. Wyoming’s data is shown in gray, as its observations are excluded due to an insufficient number of pre-treatment periods.

The missing data matrix can also be accessed via out_ub$obs.missing, where values 0, 1, 2, 3, and 4 represent missing values, control observations, treated (pre-treatment) observations, treated (post-treatment) observations, and treated observations removed from the sample, respectively.

plot(out_ub, type = "status",
     xticklabels=c(1920, 1928, 1936, 1944, 1952, 1960, 
                   1968, 1976, 1984, 1992, 2000, 2008),
     xangle=10)

Like in panelView, we can plot specific observations by setting id and adjust labels on x-axis as needed.

plot(out_ub, type = "status", xlab = "Year", ylab = "State",
     main = "Treatment Status", id = out_ub$id[out_ub$tr],
     xlim = c(1920,2012), 
     xticklabels=c(1920, 1928, 1936, 1944, 1952, 1960,
                   1968, 1976, 1984, 1992, 2000, 2008))

We re-produce the gap plot with the unbalanced panel, here, we set the range of y in between $(-10,15)$.

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

The same gap plot under the parametric bootstrap on the same unbalanced panel:

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

6.4 Bounded factor loadings

In the standard Gsynth setup, the treated-unit factor loading $\hat\lambda_i$ is obtained from an unconstrained least-squares projection of the pre-treatment residual onto the control-derived factor matrix. When few treated units lie inside the convex hull of control loadings, this projection can extrapolate aggressively, producing implied counterfactual paths that no individual control (or weighted average of them) could plausibly trace. Bounded factor loadings address this by constraining each treated loading to be a non-negative convex combination of control loadings, recovering the synthetic-control intuition that the imputed counterfactual should sit inside the support of the donor pool.

6.4.1 The simplex projection

For each treated unit $i$, the bounded estimator solves

\[ \min_{w_i \in \Delta^{N_{co}}} \;\; \frac{1}{\gamma} \, \mathrm{KL}(w_i \,\|\, \tilde{w}) \; + \; \big\| u_i^{\mathrm{pre}} - F^{\mathrm{pre}} \Lambda_{co}^\top w_i \big\|_2^2 \]

where:

$w_i$ is the $N_{co}$-vector of simplex weights for treated unit $i$, with $w_{ij} \geq 0$ and $\sum_j w_{ij} = 1$;
$\tilde{w} = (1/N_{co}, \ldots, 1/N_{co})$ is the uniform reference and $\mathrm{KL}(w \,\|\, \tilde{w}) = \sum_j w_j \log(N_{co}\, w_j)$ is the entropy regularizer;
$\gamma > 0$ controls regularization strength: as $\gamma \to \infty$ the KL term vanishes and the solution approaches a pure simplex-constrained least-squares projection; as $\gamma \to 0$ weights are forced toward uniform;
$u_i^{\mathrm{pre}}$ is the pre-treatment residual after time and unit fixed effects for unit $i$;
$F^{\mathrm{pre}}$ is the $T_0 \times r$ matrix of estimated factors restricted to the pre-treatment window;
$\Lambda_{co}$ is the $N_{co} \times r$ matrix of estimated control loadings.

The implied treated loading is $\hat\lambda_i = \Lambda_{co}^\top w_i$. Because $w_i$ lies in $\Delta^{N_{co}}$ and $\hat\lambda_i$ is therefore a convex combination of rows of $\Lambda_{co}$, the imputed counterfactual

\[\hat Y_i(0) \;=\; F \, \hat\lambda_i\]

lies pointwise in the convex hull of factor-implied control outcomes at every period.

The objective is the standard quadratic-residual fit on the pre-period plus a KL penalty toward uniform. The KL term plays two roles: it stabilizes the optimization (the unconstrained simplex projection is degenerate when $N_{co}$ exceeds the number of pre-period observations) and it shrinks toward an “equal-weight donor pool” prior, which is sensible when controls are exchangeable. The solver uses a softmax reparameterization with optim(method = "L-BFGS-B") and an analytic gradient, with a mirror-descent fallback for ill-conditioned inputs.

6.4.2 Bounded vs. unbounded: a worked comparison

The new argument loading.bound = "simplex" (default "none") activates the constraint. Regularization strength $\gamma$ is selected by 5-fold cross-validation over a log-grid by default; supply gamma.loading directly to skip CV, or pass a custom grid via gamma.loading.grid. Bounded loadings is supported for method = "gsynth" (or equivalently method = "ife", time.component.from = "nevertreated").

We fit two specifications on sim_gsynth — one unbounded (loading.bound = "none", the default) and one bounded (loading.bound = "simplex") — everything else equal:

# Unbounded gsynth (default v2.2.x behavior --- standard Xu 2017 estimator)
out.unbounded <- fect(Y ~ D, data = sim_gsynth, index = c("id", "time"),
                      method = "gsynth", r = 2,
                      se = TRUE, vartype = "parametric", nboots = 100,
                      CV = FALSE, parallel = FALSE)

# Bounded gsynth (new in v2.3.0 --- simplex projection)
out.bounded   <- fect(Y ~ D, data = sim_gsynth, index = c("id", "time"),
                      method = "gsynth", r = 2,
                      loading.bound = "simplex",
                      se = TRUE, vartype = "parametric", nboots = 100,
                      CV = FALSE, parallel = FALSE)

ATT comparison:

cat("Unbounded ATT =", round(out.unbounded$att.avg, 3), "\n")
#> Unbounded ATT = 4.64
cat("Bounded   ATT =", round(out.bounded$att.avg,   3), "\n")
#> Bounded   ATT = 5.098
cat("\ngamma.loading (CV-selected) =", round(out.bounded$gamma.loading, 4), "\n")
#> 
#> gamma.loading (CV-selected) = 100
cat("\nPer-treated projection residuals (loading.proj.resid):\n")
#> 
#> Per-treated projection residuals (loading.proj.resid):
print(round(out.bounded$loading.proj.resid, 3))
#> [1] 15.444 29.180 20.274 18.724 24.152

Large entries in loading.proj.resid (substantially above the control-fit RMSE) flag treated units near or outside the convex hull of control loadings — the simplex constraint is binding for those units, and reported intervals may under-cover (the Andrews 1999, 2001 non-standard-limit regime for constrained estimators).

6.4.2.1 Side-by-side gap plots

library(gridExtra)
gap.un  <- plot(out.unbounded, type = "gap", main = "Unbounded gsynth")
gap.bd  <- plot(out.bounded,   type = "gap", main = "Bounded gsynth (simplex)")
grid.arrange(gap.un, gap.bd, ncol = 2)

6.4.2.2 Side-by-side `loading.overlap` plots

The companion plot type loading.overlap displays factor 1 versus factor 2, with the convex hull of control loadings shaded and treated units overlaid as red triangles. The unbounded fit shows treated loadings as estimated (which can lie outside the hull); the bounded fit shows treated loadings as projected onto the simplex (always inside the hull or on its boundary).

ov.un <- plot(out.unbounded, type = "loading.overlap",
              main = "Unbounded loadings")
ov.bd <- plot(out.bounded,   type = "loading.overlap",
              main = "Bounded loadings (projected)")
grid.arrange(ov.un, ov.bd, ncol = 2)

For one-factor models ($r = 1$), the same loading.overlap diagnostic renders as a mirror histogram — treated counts above the axis, control counts below, with a vertical band marking the range of control loadings. The 1D analog of “outside the hull” is “outside the band.”

6.4.2.3 Implicit donor weights: bounded vs. unbounded

Each treated counterfactual $\hat Y_i(0) = F\,\hat\lambda_i$ is a linear combination of control-implied outcomes. The coefficients on the controls — stored in wgt.implied, an $N_{tr} \times N_{co}$ matrix — are the implicit donor weights. Under loading.bound = "none", these come from the Moore-Penrose pseudo-inverse and may be negative, exceed one, or vary wildly with small perturbations. Under loading.bound = "simplex", each row is non-negative and sums to one (a proper synthetic-control weighting). The scatter below shows every (treated unit, control unit) pair: x-axis is the unbounded weight, y-axis is the bounded weight.

library(ggplot2)
#> Warning: package 'ggplot2' was built under R version 4.5.2
W.un <- as.matrix(out.unbounded$wgt.implied)
W.bd <- as.matrix(out.bounded$wgt.implied)
df.w <- data.frame(
  unbounded = as.numeric(W.un),
  bounded   = as.numeric(W.bd)
)

ggplot(df.w, aes(x = unbounded, y = bounded)) +
  geom_hline(yintercept = 0, color = "grey80", linewidth = 0.4) +
  geom_vline(xintercept = 0, color = "grey80", linewidth = 0.4) +
  geom_abline(slope = 1, intercept = 0,
              color = "grey60", linetype = "dashed", linewidth = 0.4) +
  geom_point(alpha = 0.4, size = 1.4, color = "#3F6A99") +
  labs(x = "Unbounded weight (Moore-Penrose pseudo-inverse)",
       y = "Bounded weight (simplex-projected)",
       title = "Implicit donor weights: bounded vs unbounded",
       subtitle = "One point per (treated, control) pair; dashed line is y = x") +
  theme_bw(base_size = 11) +
  theme(panel.grid.minor = element_blank())

Three patterns are worth checking:

Many weights near $y = 0$. Bounded weights are sparse: most (treated, control) pairs receive zero weight, and only a small set of donors receives meaningful weight for each treated unit. This sparsity comes from the simplex constraint, not from the regularizer. Non-negative least squares on a probability simplex has an active-set property: at the optimum, many constraints $w_j = 0$ are active, and only a few donors have positive weight. The entropy regularizer pushes weights toward uniform $(1/N_{co})$, so it works against sparsity. With moderate $\gamma$, the constraint dominates; as $\gamma \to 0$, weights converge to uniform. Unbounded weights, by contrast, spread non-zero mass across the donor pool because the pseudo-inverse has no non-negativity constraint and no preference for sparsity.
Negative unbounded weights become zero. Points with $x < 0$ sit at $y = 0$, or close to it. The simplex projection cannot use negative weights, so it removes those donors from the counterfactual. This matches the synthetic-control intuition: a counterfactual should not be a negative combination of donors.
Bounded weights cannot exceed one. Points with $x > 1$ move to $y \leq 1$. When the pseudo-inverse puts too much weight on one donor, the simplex caps that weight and redistributes the rest across other donors.

The result is a different counterfactual representation. The pseudo-inverse is dense and can extrapolate. The simplex is sparse, uses interpolation only, and gives each weight a synthetic-control interpretation.

6.4.2.4 Reading the comparison

Because sim_gsynth follows a clean IFE DGP and treated units sit reliably inside the convex hull of control loadings, the two fits give very similar ATTs — the simplex constraint barely binds at the loading level. The weight representations differ substantially even when the implied counterfactuals do not, which is informative on its own: bounded loadings buy interpretability and robustness without paying much in fit on well-behaved data. Real datasets with sparse donor pools (few treated units alongside controls with disparate factor loadings) typically show larger ATT gaps; in those settings, the simplex bound corrects the extrapolation that pure least-squares performs and provides a more honest synthetic-control counterfactual.

6.4.3 Caveats

Currently restricted to method = "ife" or "gsynth" (equivalently, time.component.from = "nevertreated").
Percentile-bootstrap intervals may under-cover when the simplex constraint binds (true treated loading on the boundary of $\mathrm{conv}(\Lambda_{co})$). Detect via loading.proj.resid. Boundary-corrected inference is deferred to a later release.
The unit fixed-effect scalar $\alpha_i$ is not bounded; under force %in% c("unit", "two-way") it is computed as the residual mean per treated unit.

6.5 CFE with GSynth

Gsynth is not limited to the IFE method. The complex fixed effects (CFE) estimator can also be used with time.component.from = "nevertreated". This is useful when the data generating process involves not only latent interactive fixed effects but also additional additive fixed effects (e.g., region, industry), unit-specific time trends (via Q.type), or time-invariant covariates with time-varying coefficients (Z/gamma). See Chapter 5 for the full CFE tutorial.

6.5.1 Example: CFE + nevertreated

out.cfe.nt <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id","time"),
                   method = "cfe", force = "two-way",
                   time.component.from = "nevertreated",
                   Q.type = "linear",
                   se = FALSE, CV = TRUE, r = c(0, 5),
                   parallel = TRUE, cores = 16)

cat("CFE + nevertreated: r.cv =", out.cfe.nt$r.cv,
    ", ATT =", round(out.cfe.nt$att.avg, 3), "\n")
#> CFE + nevertreated: r.cv = 2 , ATT = 6.024

Under time.component.from = "nevertreated", the CFE estimator uses only never-treated controls for factor estimation — the same principle as gsynth. The additional CFE components (extra FE, time trends, Z/gamma) are estimated on the control panel and projected onto treated units alongside the interactive factors. See Chapter 5 for the full CFE tutorial.

6.5.2 Comparing gsynth and CFE + nevertreated

Both approaches share the same estimation regime (never-treated controls only), but CFE can add structural components on top. Whether the extra flexibility helps depends on the data generating process. We compare three specifications using fect_mspe():

# Model 1: gsynth (pure IFE, r = 2)
out.gsynth.comp <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id","time"),
                        method = "gsynth", force = "two-way",
                        r = 2, se = FALSE, CV = FALSE, max.iteration = 20000)

# Model 2: CFE + nevertreated with r = 2 only (equivalent to gsynth)
out.cfe.nt.comp <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id","time"),
                        method = "cfe", force = "two-way",
                        time.component.from = "nevertreated",
                        r = 2, se = FALSE, CV = FALSE, max.iteration = 20000)

# Model 3: CFE + nevertreated with r = 2 and linear trend (overspecified)
out.cfe.nt.lin <- fect(Y ~ D + X1 + X2, data = sim_gsynth, index = c("id","time"),
                       method = "cfe", force = "two-way",
                       time.component.from = "nevertreated",
                       Q.type = "linear", r = 2, se = FALSE, CV = FALSE,
                       max.iteration = 20000)

mspe.comp <- fect_mspe(list(gsynth_r2 = out.gsynth.comp,
                            CFE_r2 = out.cfe.nt.comp,
                            CFE_linear_r2 = out.cfe.nt.lin), seed = 1234)
print(mspe.comp$summary[, c("Model", "MSPE", "RMSE", "MAD")])
#>           Model     MSPE     RMSE      MAD
#> 1 CFE_linear_r2 110.3310 10.50386 58.19996
#> 2        CFE_r2 110.3310 10.50386 58.19996
#> 3     gsynth_r2 100.9605 10.04791 45.76502

Since sim_gsynth follows a pure IFE data generating process with two factors, Models 1 and 2 should produce identical MSPE — confirming that CFE with time.component.from = "nevertreated" and no additional structure is numerically equivalent to gsynth. Model 3, which adds unnecessary linear trends, should produce similar or slightly worse MSPE because the extra parameters add noise without benefit when the true DGP has no unit-specific trends.

6.6 Inference

The default vartype for the gsynth-style regime is "parametric" — a two-stage pseudo-treated bootstrap that targets the conditional variance $V_t = \mathrm{Var}(\widehat{\mathrm{ATT}}_t - \mathrm{ATT}_t \mid \Lambda, F, X, D)$ on small-$N_{\text{tr}}$ panels where the standard nonparametric bootstrap is unreliable. The para.error argument (introduced in v2.4.2) chooses how the residual error process is resampled: "auto" (default, resolves to "empirical" on a fully-observed panel and "ar" on a missing-data panel), "ar" (the v2.4.1 behavior), "empirical", or "wild". The (vartype × ci.method × para.error) interaction, parametric-bootstrap mechanics, and an empirical coverage study live in Chapter 7, which is the natural next read after this chapter.

6.7 Additional Notes

Unbalanced Panels: Running Gsynth within fect on unbalanced panels will take significantly more time compared to balanced panels, often by a factor of 100:1 or more. This is because the EM algorithm, which fills in missing values (implemented in C++), requires many more iterations to converge. To reduce run-time, users can remove units or time periods with extensive missing values. Understanding the data structure before running any regressions is always helpful. Note that this approach differs from setting method = "ife", as no pre-treatment data from the treated units is used to impute $\hat{Y}(0)$.
Adding Covariates: Including covariates in the model will significantly slow down the algorithm, as the IFE/MC model requires more time to converge. Users should be aware of this trade-off when incorporating covariates.
Setting min.T0: Setting min.T0 to a positive value helps. The algorithm will automatically exclude treated units with too few pre-treatment periods. A larger $T_0$ reduces bias in causal estimates and minimizes the risk of severe extrapolation. When running cross-validation to select the number of factors, min.T0 must be equal to or greater than (r.max + 2). Errors frequently occur when there are too few pre-treatment periods, so ensuring adequate $T_0$ (e.g. setting min.T0 = 5) is crucial.

6.8 How to Cite

If you find this method helpful, you can cite Xu (2017).

@article{Xu2017,
  title = {Generalized Synthetic Control Method: Causal Inference with Interactive Fixed Effects Models},
  author = {Xu, Yiqing},
  journal = {Political Analysis},
  volume = {25},
  number = {1},
  pages = {57--76},
  year = {2017}
}

6.1 Basic Usage

6.1.1 Simulated Data

6.1.2 Estimation

6.1.3 Cross-validation

6.1.4 Rolling-window CV

6.1.5 Visualizing Results

6.2 Staggered DID

6.2.1 Data structure

6.2.2 Estimation w/o factors

6.2.3 Estimation w/ factors

6.2.4 Implied weights

6.2.5 Visualizing results

6.2.6 Cumulative Effects

6.3 Unbalanced panels

6.3.1 Estimation

6.3.2 Parametric bootstrap on an unbalanced panel

6.3.3 Visualizing results

6.4 Bounded factor loadings

6.4.1 The simplex projection

6.4.2 Bounded vs. unbounded: a worked comparison

6.4.2.1 Side-by-side gap plots

6.4.2.2 Side-by-side loading.overlap plots

6.4.2.3 Implicit donor weights: bounded vs. unbounded

6.4.2.4 Reading the comparison

6.4.3 Caveats

6.5 CFE with GSynth

6.5.1 Example: CFE + nevertreated

6.5.2 Comparing gsynth and CFE + nevertreated

6.6 Inference

6.7 Additional Notes

6.8 How to Cite

6.4.2.2 Side-by-side `loading.overlap` plots