(Some) modern staggered difference-in-differences estimators for Julia.
StagDiDModels.jl implements some of the recent advances in difference-in-differences estimation that address heterogeneous treatment effects in staggered treatment adoption. All estimators provide standard errors, model diagnostics, and are compatible with the StatsAPI ecosystem. Integrates nicely with HonestDiD.jl
| Estimator | Function | Reference |
|---|---|---|
| BJS Imputation | fit_bjs() |
Borusyak, Jaravel, Spiess (2023) |
| Gardner Two-Stage | fit_gardner_static(), fit_gardner_dynamic() |
Gardner (2021) |
| Sun-Abraham | fit_sunab() |
Sun & Abraham (2021) |
| TWFE | fit_twfe_static(), fit_twfe_dynamic() |
Standard approach |
using Pkg
Pkg.add(url="https://github.com/eohne/StagDiDModels.jl")
using StagDiDModels
using DataFramesPanel data with the following variables:
- Unit identifier (e.g., firm, individual)
- Time period (integer, e.g., year, quarter)
- Treatment timing (when unit first receives treatment)
- Outcome variable
Never-treated units must be coded as g = 0 or g = missing.
Static average treatment effect:
julia> fit_bjs_static(df; y=:dep_var, id=:unit, t=:year, g=:g,cluster=:unit)
BJSModel
==========================================
Number of obs: 31000 Outcome: :dep_var
Treated obs: 10253 Est. Type: :static
Donor obs: 20747 R² First Stage:0.777
Cluster: :unit
==========================================
Coef. Std. Error t Pr(>|t|)
──────────────────────────────────────────
_ATT 2.26295 0.0313968 72.08 <1e-99
==========================================Dynamic event study:
fit_bjs_dynamic(df; y=:dep_var, id=:unit, t=:year, g=:g,
horizons=true, pretrends=true, cluster=:unit)
┌ Warning: Collinear pre-trends omitted: pre20
└ @ Main c:\Users\eohneberg\OneDrive\Julia\Packages\DiDModels.jl\src\estimators\bjs.jl:492
BJSModel
================================================
Number of obs: 31000 Outcome: :dep_var
Treated obs: 10253 Est. Type: :dynamic
Donor obs: 20747 R² First Stage: 0.777
Cluster: :unit
================================================
Coef. Std. Error t Pr(>|t|)
────────────────────────────────────────────────
τ::-19 0.00494698 0.142763 0.03 0.9724
τ::-18 -0.094598 0.142099 -0.67 0.5057
τ::-17 -0.0812596 0.144141 -0.56 0.5731
τ::-16 -0.112528 0.13729 -0.82 0.4126
τ::-15 -0.0753304 0.14184 -0.53 0.5955
τ::-14 0.05654 0.144027 0.39 0.6947
τ::-13 -0.068671 0.143188 -0.48 0.6316
τ::-12 -0.0915126 0.14232 -0.64 0.5204
τ::-11 -0.243579 0.142098 -1.71 0.0868
τ::-10 -0.077792 0.106634 -0.73 0.4659
τ::-9 -0.052073 0.129615 -0.40 0.6880
τ::-8 -0.123991 0.128756 -0.96 0.3358
τ::-7 -0.00807412 0.129888 -0.06 0.9504
τ::-6 -0.242521 0.131568 -1.84 0.0656
τ::-5 -0.157734 0.130259 -1.21 0.2262
τ::-4 -0.0704127 0.126895 -0.55 0.5791
τ::-3 -0.105707 0.132233 -0.80 0.4242
τ::-2 -0.0469036 0.13182 -0.36 0.7221
τ::-1 -0.11832 0.128169 -0.92 0.3561
τ::0 1.51314 0.0754774 20.05 <1e-74
τ::1 1.66384 0.0767514 21.68 <1e-84
τ::2 1.86437 0.0745015 25.02 <1e-99
τ::3 1.91872 0.074717 25.68 <1e-99
τ::4 1.87322 0.0741817 25.25 <1e-99
τ::5 1.87845 0.0756719 24.82 <1e-99
τ::6 2.14373 0.0763269 28.09 <1e-99
τ::7 2.23778 0.0761084 29.40 <1e-99
τ::8 2.3365 0.0744627 31.38 <1e-99
τ::9 2.34353 0.0747168 31.37 <1e-99
τ::10 2.53443 0.0810955 31.25 <1e-99
τ::11 2.47945 0.119535 20.74 <1e-79
τ::12 2.63494 0.115318 22.85 <1e-92
τ::13 2.9445 0.110473 26.65 <1e-99
τ::14 2.78171 0.114664 24.26 <1e-99
τ::15 2.71471 0.120305 22.57 <1e-90
τ::16 2.88065 0.115632 24.91 <1e-99
τ::17 2.99384 0.114385 26.17 <1e-99
τ::18 2.64617 0.115458 22.92 <1e-93
τ::19 2.87531 0.114058 25.21 <1e-99
τ::20 2.90466 0.113202 25.66 <1e-99
================================================General specification:
julia> model = fit_bjs(df;
y=:dep_var, # Outcome variable
id=:unit, # Unit identifier
t=:year, # Time variable
g=:g, # Treatment timing
controls=Symbol[], # Control variables
fe=(:unit, :year), # Fixed effects (default: unit + time)
weights=:w, # Observation weights
cluster=:unit, # Clustering variable
horizons=[0,1,2,3,4], # Post-treatment periods
pretrends=[-3,-2,-1], # Pre-treatment periods
)
BJSModel
===============================================
Number of obs: 31000 Outcome: :dep_var
Treated obs: 10253 Est. Type: :dynamic
Donor obs: 20747 R² First Stage: 0.780
Cluster: :unit
===============================================
Coef. Std. Error t Pr(>|t|)
───────────────────────────────────────────────
τ::-3 -0.00691821 0.0829011 -0.08 0.9335
τ::-2 0.0415102 0.0828488 0.50 0.6165
τ::-1 0.0295101 0.0802191 0.37 0.7130
τ::0 1.53155 0.0828972 18.48 <1e-65
τ::1 1.61755 0.0845701 19.13 <1e-69
τ::2 1.86139 0.0822017 22.64 <1e-91
τ::3 1.97443 0.0811454 24.33 <1e-99
τ::4 1.82614 0.079427 22.99 <1e-93
===============================================Sometimes you expect the treatment effect to differ along some dimension. There
are two ways to look at this, both matching the Stata did_imputation options of
the same name.
Use project when the dimension is continuous and you want to know how the effect
changes with it. It regresses the individual imputed effects on the covariates you
pass and reports an intercept (τ_cons) together with a slope for each covariate
(τ_<var>). A positive, significant slope means the effect is larger when that
variable is larger. With horizons you get one intercept and slope per event time.
# Does the effect grow with x1?
fit_bjs(df; y=:dep_var, id=:unit, t=:year, g=:g, project=[:x1], cluster=:unit)Use hetby when the dimension is discrete (up to 30 groups). It reports a separate
effect for each group value, named τ_<value>, or τ<h>_<value> when combined with
horizons.
# A separate ATT for each region
fit_bjs(df; y=:dep_var, id=:unit, t=:year, g=:g, hetby=:region, cluster=:unit)project and hetby cannot be used together.
The standard errors on project slopes and hetby group effects are not naive
OLS errors on the imputed effects. They use the same clustered influence function
variance that the package applies to every other treatment coefficient, and they
honor the cluster argument. So for project the slope inference is clustered at
the level you pass, and you do not need to cluster or bootstrap it by hand.
The usual caveat applies: the errors are only as reliable as the number of clusters,
so be careful with very few clusters.
minn sets the minimum effective sample size a coefficient needs before it is
reported. It defaults to 0, which keeps every coefficient, so existing code is
unaffected. Setting minn=30 reproduces Stata's default behaviour and drops
coefficients (often long horizons or small heterogeneity cells) that rest on too
few effective observations, since their standard errors can be misleadingly small.
Suppressed coefficients are reported in a warning.
fit_bjs(df; y=:dep_var, id=:unit, t=:year, g=:g, horizons=true, minn=30, cluster=:unit)Static:
julia> fit_gardner_static(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit)
GardnerModel
==========================================
Number of obs: 31000 R²: 0.353
Treated obs: 10253 R² adjusted: 0.353
Donor obs: 20747 F-statistic: 16904.0
Outcome: :dep_var Cluster : :unit
==========================================
Coef. Std. Error t Pr(>|t|)
──────────────────────────────────────────
_ATT 2.26295 0.0338789 66.80 <1e-99
==========================================Dynamic:
julia> fit_gardner_dynamic(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit, ref_p=-1)
GardnerModel
=================================================
Number of obs: 31000 R²: 0.373
Treated obs: 10253 R² adjusted: 0.372
Donor obs: 20747 F-statistic: 460.171
Outcome: :dep_var Cluster : :unit
=================================================
Coef. Std. Error t Pr(>|t|)
─────────────────────────────────────────────────
τ::-20 0.0482675 0.0604799 0.80 0.4248
τ::-19 0.0434281 0.0607702 0.71 0.4748
τ::-18 -0.000894858 0.0609516 -0.01 0.9883
τ::-17 -0.028753 0.0632599 -0.45 0.6495
τ::-16 0.0246948 0.0602951 0.41 0.6821
τ::-15 0.0228411 0.0612728 0.37 0.7093
τ::-14 0.0840854 0.064303 1.31 0.1910
τ::-13 0.0107916 0.0615972 0.18 0.8609
τ::-12 -0.023326 0.0637062 -0.37 0.7143
τ::-11 -0.10326 0.0639668 -1.61 0.1065
τ::-10 -0.00170496 0.0413833 -0.04 0.9671
τ::-9 0.0129794 0.0411607 0.32 0.7525
τ::-8 -0.0143057 0.0416646 -0.34 0.7313
τ::-7 0.0533169 0.0409899 1.30 0.1934
τ::-6 -0.082763 0.0412686 -2.01 0.0449
τ::-5 -0.0377117 0.0413368 -0.91 0.3616
τ::-4 -0.00635132 0.0407193 -0.16 0.8761
τ::-3 -0.00746851 0.0414751 -0.18 0.8571
τ::-2 0.0315947 0.0421007 0.75 0.4530
τ::0 1.51314 0.076538 19.77 <1e-85
τ::1 1.66384 0.0782553 21.26 <1e-98
τ::2 1.86437 0.0761971 24.47 <1e-99
τ::3 1.91872 0.0759098 25.28 <1e-99
τ::4 1.87322 0.0750175 24.97 <1e-99
τ::5 1.87845 0.0765166 24.55 <1e-99
τ::6 2.14373 0.0763666 28.07 <1e-99
τ::7 2.23778 0.0763247 29.32 <1e-99
τ::8 2.3365 0.0749124 31.19 <1e-99
τ::9 2.34353 0.0747161 31.37 <1e-99
τ::10 2.53443 0.0811948 31.21 <1e-99
τ::11 2.47945 0.119535 20.74 <1e-94
τ::12 2.63494 0.115318 22.85 <1e-99
τ::13 2.9445 0.110473 26.65 <1e-99
τ::14 2.78171 0.114664 24.26 <1e-99
τ::15 2.71471 0.120305 22.57 <1e-99
τ::16 2.88065 0.115632 24.91 <1e-99
τ::17 2.99384 0.114385 26.17 <1e-99
τ::18 2.64617 0.115458 22.92 <1e-99
τ::19 2.87531 0.114058 25.21 <1e-99
τ::20 2.90466 0.113202 25.66 <1e-99
=================================================ATT:
julia> fit_sunab(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit, agg=:att)
SunabModel
=========================================
Number of obs:31000 R²: 0.799
Treated obs: 10253 R² adjusted: 0.792
Pre obs: 20747 F-statistic:1.23502e5
Dep. Var: :dep_var Cluster: :unit
=========================================
Coef. Std. Error t Pr(>|t|)
─────────────────────────────────────────
ATT 2.32119 0.0820088 28.30 <1e-99
=========================================Dynamic coefficients:
julia> fit_sunab(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit, agg=:dynamic)
SunabModel
=================================================
Number of obs: 31000 R²: 0.799
Treated obs: 10253 R² adjusted: 0.792
Pre obs: 20747 F-statistic: 3083.66
Dep. Var: :dep_var Cluster: :unit
=================================================
Coef. Std. Error t Pr(>|t|)
─────────────────────────────────────────────────
τ::-20 -0.0151877 0.147769 -0.10 0.9182
τ::-19 -0.0082711 0.147491 -0.06 0.9553
τ::-18 -0.0997613 0.152664 -0.65 0.5136
τ::-17 -0.0760892 0.151675 -0.50 0.6160
τ::-16 -0.113553 0.152684 -0.74 0.4572
τ::-15 -0.121286 0.152397 -0.80 0.4263
τ::-14 0.000968567 0.150746 0.01 0.9949
τ::-13 -0.0875309 0.148901 -0.59 0.5568
τ::-12 -0.0281754 0.154189 -0.18 0.8550
τ::-11 -0.324711 0.152768 -2.13 0.0338
τ::-10 0.0398103 0.108042 0.37 0.7126
τ::-9 0.0655079 0.108619 0.60 0.5466
τ::-8 -0.00649778 0.107914 -0.06 0.9520
τ::-7 0.109307 0.108343 1.01 0.3133
τ::-6 -0.125073 0.108768 -1.15 0.2505
τ::-5 -0.039797 0.108363 -0.37 0.7135
τ::-4 0.0476292 0.10436 0.46 0.6482
τ::-3 0.0119353 0.103833 0.11 0.9085
τ::-2 0.0698439 0.109226 0.64 0.5227
τ::0 1.53545 0.109545 14.02 <1e-40
τ::1 1.69124 0.109982 15.38 <1e-47
τ::2 1.87083 0.106399 17.58 <1e-59
τ::3 1.9474 0.108805 17.90 <1e-61
τ::4 1.84971 0.106231 17.41 <1e-58
τ::5 1.89713 0.109857 17.27 <1e-57
τ::6 2.18776 0.107846 20.29 <1e-76
τ::7 2.25533 0.105809 21.32 <1e-82
τ::8 2.32726 0.108476 21.45 <1e-83
τ::9 2.3889 0.0952692 25.08 <1e-99
τ::10 2.55571 0.109845 23.27 <1e-95
τ::11 2.62696 0.161015 16.32 <1e-52
τ::12 2.78245 0.153434 18.13 <1e-63
τ::13 3.09201 0.153139 20.19 <1e-75
τ::14 2.92923 0.152714 19.18 <1e-69
τ::15 2.86222 0.158084 18.11 <1e-62
τ::16 3.02817 0.155266 19.50 <1e-71
τ::17 3.14135 0.152887 20.55 <1e-77
τ::18 2.79368 0.158095 17.67 <1e-60
τ::19 3.02282 0.153328 19.71 <1e-72
τ::20 3.05217 0.152721 19.99 <1e-74
=================================================Static:
julia> fit_twfe_static(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit)
TWFEModel
==========================================
Number of obs: 31000 R²: 0.794
Post obs: 21320 R² adjusted: 0.787
Pre obs: 9680 F-statistic:119798.0
Dep. Var: :dep_var Cluster: :unit
==========================================
Coef. Std. Error t Pr(>|t|)
──────────────────────────────────────────
post 2.01215 0.0311817 64.53 <1e-99
==========================================
Dynamic:
julia> fit_twfe_dynamic(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit, ref_p=-1)
TWFEModel
================================================
Number of obs: 31000 R²: 0.799
Post obs: 10253 R² adjusted: 0.792
Pre obs: 20747 F-statistic: 3073.42
Dep. Var: :dep_var Cluster: :unit
================================================
Coef. Std. Error t Pr(>|t|)
────────────────────────────────────────────────
τ::-20 0.228396 0.121549 1.88 0.0605
τ::-19 0.225683 0.122522 1.84 0.0658
τ::-18 0.124711 0.127555 0.98 0.3285
τ::-17 0.143929 0.127102 1.13 0.2577
τ::-16 0.122706 0.123789 0.99 0.3218
τ::-15 0.159711 0.126671 1.26 0.2077
τ::-14 0.317349 0.126036 2.52 0.0120
τ::-13 0.181699 0.123417 1.47 0.1413
τ::-12 0.148058 0.127864 1.16 0.2472
τ::-11 0.0253997 0.104909 0.24 0.8087
τ::-10 -0.0456146 0.104611 -0.44 0.6629
τ::-9 -0.0361548 0.106546 -0.34 0.7344
τ::-8 -0.111102 0.105043 -1.06 0.2905
τ::-7 0.0172967 0.105776 0.16 0.8701
τ::-6 -0.195828 0.105384 -1.86 0.0634
τ::-5 -0.11145 0.105431 -1.06 0.2907
τ::-4 0.0305665 0.101692 0.30 0.7638
τ::-3 -0.0268863 0.10111 -0.27 0.7904
τ::-2 0.00899443 0.106906 0.08 0.9330
τ::0 1.47408 0.104917 14.05 <1e-40
τ::1 1.63928 0.107153 15.30 <1e-46
τ::2 1.82363 0.103913 17.55 <1e-59
τ::3 1.90212 0.106412 17.87 <1e-61
τ::4 1.79594 0.10399 17.27 <1e-57
τ::5 1.82125 0.107351 16.97 <1e-56
τ::6 2.09319 0.104439 20.04 <1e-74
τ::7 2.18456 0.103711 21.06 <1e-81
τ::8 2.30298 0.106062 21.71 <1e-85
τ::9 2.27728 0.091008 25.02 <1e-99
τ::10 2.55733 0.108288 23.62 <1e-97
τ::11 2.44798 0.136409 17.95 <1e-61
τ::12 2.6095 0.131919 19.78 <1e-73
τ::13 2.89423 0.130662 22.15 <1e-88
τ::14 2.68905 0.127594 21.08 <1e-81
τ::15 2.62285 0.135489 19.36 <1e-70
τ::16 2.68003 0.132898 20.17 <1e-75
τ::17 2.83728 0.128553 22.07 <1e-87
τ::18 2.5352 0.133589 18.98 <1e-68
τ::19 2.6402 0.125345 21.06 <1e-81
τ::20 2.76279 0.129787 21.29 <1e-82
================================================All models support standard StatsAPI functions:
coef(model) # Coefficient estimates
vcov(model) # Variance-covariance matrix
stderror(model) # Standard errors
coefnames(model) # Coefficient names
confint(model) # Confidence intervals
nobs(model) # Number of observationsjulia> model = fit_sunab(df; y=:dep_var, id=:unit, t=:year, g=:g,
cluster=:unit, agg=:dynamic);
julia> cum_effects = cumulative_effects(model);
# Access results
julia> cum_effects.τ # Event times
40-element Vector{Int64}:
-20
-19
⋮
19
20
julia> cum_effects.cumulative # Cumulative effects
40-element Vector{Float64}:
-0.7009304415811053
-0.6857427528313305
⋮
48.78563490491884
51.83780653599853
julia> cum_effects.std_errors # Standard errors
40-element Vector{Float64}:
1.6669189569635976
1.5663756814986112
⋮
1.728900903879408
1.8318442896860898using Makie, CairoMakie
m = fit_bjs_dynamic(df; y=:dep_var, id=:unit, t=:year, g=:g,cluster=:unit);
fig = plot_event_study(m; title="BJS Event Study", color=:steelblue)result = plot_comparison(df;
y=:dep_var, id=:unit, t=:year, g=:g,
models=[:bjs, :gardner, :sunab, :twfe],
cluster=:unit).figureAll estimator functions accept:
df::DataFrame: Panel datasety::Symbol: Outcome variable nameid::Symbol: Unit identifiert::Symbol: Time period variableg::Symbol: Treatment timing variablecontrols::Vector{Symbol}: Additional controls (default:Symbol[])cluster::Symbol: Clustering variable (default: clusters byid)weights::Union{Nothing,Symbol}: Observation weights (default:nothing)
fit_bjs additionally accepts:
project::Vector{Symbol}: Continuous treatment effect heterogeneity, the Stataproject()option (default:Symbol[])hetby::Union{Nothing,Symbol}: Discrete treatment effect heterogeneity, the Statahetby()option (default:nothing)minn::Real: Minimum effective sample size per coefficient.0keeps everything,30matches Stata's default (default:0)
| StagDiDModels.jl | R | Stata |
|---|---|---|
fit_bjs, fit_bjs_dynamic, fit_bjs_static |
didimputation::did_imputation † |
did_imputation |
fit_gardner_dynamic, fit_gardner_static |
did2s::did2s |
— |
fit_sunab |
fixest::feols(y ~ sunab(g, t) | id + t) |
— |
fit_twfe_dynamic, fit_twfe_static |
fixest::feols(y ~ i(event_time) | id + t) |
— |
† The R didimputation package reports conservative standard errors (it omits the
influence-function variance) and lacks the hetby / project options;
StagDiDModels follows the Stata did_imputation implementation. See Performance.
All four estimators are pure Julia and were benchmarked against the standard reference implementations on synthetic staggered-adoption panels. Timings are warm (compilation excluded); reference tools and the BJS "single-threaded" row run on one thread for fairness, on an Intel Core Ultra 9 185H (Julia 1.11, R 4.5.1, StataNow 19 SE). Across all four estimators, coefficients and standard errors reproduce the reference implementations to ~7–8 significant figures — down to floating-point differences between regression backends.
120,000 observations (10,000 units × 12 periods), dynamic event study with clustered standard errors:
| Estimator | Reference | Reference | StagDiDModels | Speed-up |
|---|---|---|---|---|
| BJS imputation — single-threaded | Stata did_imputation |
12.5 s | 1.30 s | 9.6× |
| BJS imputation — multithreaded (16 cores) | Stata did_imputation |
12.5 s | 0.93 s | 13.4× |
| Gardner two-stage | R did2s |
18.4 s | 4.02 s | 4.6× |
| Sun–Abraham | R fixest::sunab |
0.16 s | 0.12 s | 1.3× |
| TWFE event study | R fixest::feols |
0.04 s | 0.08 s | 0.5× |
Notes:
- BJS is ~9–11× faster than Stata
did_imputationsingle-threaded and ~12–15× with the multithreaded path (multithreaded=true, the default), roughly flat across 12k–120k obs. Stata SE is single-threaded; even granting Stata MP an optimistic 2–3× the lead holds. - Gardner speed-up widens on smaller panels — ~16× at 12k obs, narrowing
to 4.6× at 120k because its influence-function SE solve has fixed-effect
fill-in that grows with the number of units. (
did2sscales worse still.) - Sun–Abraham is a modest win (~1.3× at 120k, ~2.5× at 12k); it's a single saturated regression plus a cheap aggregation, so it's backend-bound.
- TWFE is effectively a tie:
fixest's compiled core is marginally faster on a plain event-studyfeols. The value here is the unified API and matching output, not raw speed. - The R
didimputationpackage runs the BJS estimator faster than Stata, but it reports conservative standard errors (it omits the influence-function variance) and lacks thehetby/projectheterogeneity options. StagDiDModels follows the Statadid_imputationimplementation — which it reproduces exactly — sodidimputationis not a like-for-like comparison and is not shown.
Benchmark harnesses are in bench/.
- Borusyak, K., Jaravel, X., & Spiess, J. (2023). Revisiting event study designs: Robust and efficient estimation.
- Gardner, J. (2022). Two-stage differences in differences.
- Sun, L., & Abraham, S. (2021). Estimating dynamic treatment effects in event studies with heterogeneous treatment effects. Journal of Econometrics, 225(2), 175-199.