*! ═════════════════════════════════════════════════════════════════ *! Double LASSO in Stata: replicating Donohue & Levitt (2001) on *! abortion and crime. *! *! Stata companion to the R tutorial in content/tutorials/r_double_lasso/. *! Same data, same five estimators, same identification story. *! *! 1. First-difference OLS (no controls) *! 2. Kitchen-sink OLS (all 284 partialled controls) *! 3. PSL (Post-Structural LASSO) (one rlasso, treatment pinned via pnotpen) *! 4. DL-rigorous (Belloni-Chernozhukov-Hansen penalty via -rlasso-) *! 5. DL-CV (cross-validated penalty via -cvlasso-, lcount(10)) *! *! Engines: rlasso / cvlasso / lasso2 / pdslasso from the StataLasso *! suite (Ahrens, Hansen, Schaffer 2018, "pdslasso and ivlasso"). *! *! Usage: stata -b do analysis.do *! Outputs: analysis.log *! results_table2.csv (5 estimators x 3 outcomes) *! selection_diagnostic.csv (|I_y|, |I_d|, |union|) *! stata_double_lasso_*.png (3 dark-theme figures: forest plot, *! selection bars, rigorous-vs-CV compare) *! *! Data: Six CSVs pulled over HTTPS from the r_double_lasso post, *! which itself extracted them from the JAE replication *! archive's .mat files via prepare_data.R. *! *! References: *! - Fitzgerald, Lattimore, Robinson & Zhu (2026), JAE. *! - Belloni, Chernozhukov & Hansen (2014), Rev. Econ. Stud. 81: 608-650. *! - Donohue & Levitt (2001), Q. J. Econ. 116: 379-420. *! - Ahrens, Hansen & Schaffer (2018), pdslasso package. *! ═════════════════════════════════════════════════════════════════ * ── 0. Setup ───────────────────────────────────────────────────── version 18 clear all set more off set seed 20260520 capture log close log using "analysis.log", replace text * Check required packages are installed (already in user's ado/plus). foreach pkg in rlasso cvlasso lasso2 pdslasso coefplot { capture which `pkg' if _rc { di as error "Required package not installed: `pkg'" di as error "Install via: ssc install lassopack" di as error " ssc install pdslasso" di as error " ssc install coefplot" exit 198 } } * Site palette (dark theme, mirrors R version's ggplot theme_site). * DARK_BG = #0f1729 = "15 23 41" * DARK_PNL = #1f2b5e = "31 43 94" * LIGHT_TX = #c8d0e0 = "200 208 224" * LIGHTER = #e8ecf2 = "232 236 242" * STEEL = #6a9bcc = "106 155 204" * ORANGE = #d97757 = "217 119 87" * TEAL = #00d4c8 = "0 212 200" * LT_ORG = #e8956a = "232 149 106" global C_BG "15 23 41" global C_PANEL "31 43 94" global C_GRID "31 43 94" global C_TEXT "200 208 224" global C_TXTHI "232 236 242" global C_STEEL "106 155 204" global C_ORG "217 119 87" global C_TEAL "0 212 200" global C_LTORG "232 149 106" * Standard graph-region options for the dark theme. global DARKBG /// graphregion(fcolor("$C_BG") ifcolor("$C_BG") lcolor("$C_BG") ilcolor("$C_BG")) /// plotregion(fcolor("$C_BG") lcolor("$C_BG")) * ── 1. Data loading (six CSVs from GitHub raw URLs) ─────────────── di _n(2) "========================================" di "STEP 1 - DATA LOADING (six CSVs over HTTPS)" di "========================================" local BASE = "https://raw.githubusercontent.com/cmg777/starter-academic-v501/master/content/tutorials/r_double_lasso/data" tempfile linear partialled ctrl_v ctrl_p ctrl_m * (a) Raw first-differenced outcomes and treatments + state IDs. * Cols: state, Dyv, Dxv, Dyp, Dxp, Dym, Dxm import delimited "`BASE'/levitt_linear.csv", clear varnames(1) case(preserve) gen long obs_id = _n di " levitt_linear.csv : `=_N' obs x `=c(k)' cols" save "`linear'" * (b) Partialled (year-FE-removed) outcomes and treatments. * Cols: state, DxV, DyV, DxP, DyP, DxM, DyM import delimited "`BASE'/levitt_partialled.csv", clear varnames(1) case(preserve) * drop the duplicate state col so the merge is clean drop state gen long obs_id = _n di " levitt_partialled.csv : `=_N' obs x `=c(k)' cols (incl obs_id)" save "`partialled'" * (c)-(e) Three 284-column control matrices, one per outcome. * Column names in source CSV use ^, *, ( ) etc. which Stata * sanitises on import; we then rename to zv1..zv284 (etc.) * so downstream code can address them uniformly. foreach o in v p m { if "`o'" == "v" local long "viol" if "`o'" == "p" local long "prop" if "`o'" == "m" local long "murd" import delimited "`BASE'/levitt_controls_`long'.csv", clear varnames(1) local k = 0 foreach var of varlist _all { local ++k rename `var' z`o'`k' } di " levitt_controls_`long'.csv : `=_N' obs x `k' cols (renamed zv1..z`o'`k')" gen long obs_id = _n save "`ctrl_`o''" } * Combine into one working dataset. use "`linear'", clear merge 1:1 obs_id using "`partialled'", nogen merge 1:1 obs_id using "`ctrl_v'", nogen merge 1:1 obs_id using "`ctrl_p'", nogen merge 1:1 obs_id using "`ctrl_m'", nogen * Sanity check: must be 576 obs, 48 states. assert _N == 576 qui levelsof state, local(states) local nstates : word count `states' assert `nstates' == 48 di _n " Merged working dataset: `=_N' obs, `nstates' clusters (states)" * ── 2. Convenience: outcome metadata ───────────────────────────── * For each outcome we have: * raw differenced: y = Dy[v/p/m] d = Dx[v/p/m] * year-FE partialled:y = Dy[V/P/M] d = Dx[V/P/M] * 284 partialled controls: z[v/p/m]1..z[v/p/m]284 * * The partialling step (done in the original Matlab pre-processing) * absorbs year fixed effects via Frisch-Waugh-Lovell: every variable * v becomes v - T (T'T)^-1 T' v, where T is the matrix of year * dummies. Regressions on the partialled vars are equivalent to * regressing the raw differences while controlling for year dummies, * with one less degree of freedom. * Will use these foreach loops repeatedly. * Iteration uses short prefixes v/p/m (matching the renamed control vars * zv1..zv284, zp1..zp284, zm1..zm284). Look-up helpers below map the * short prefix to outcome label and raw/partialled variable names. * * prefix -> label y (raw) d (raw) y (partialled) d (partialled) * v -> Violent crime Dyv Dxv DyV DxV * p -> Property crime Dyp Dxp DyP DxP * m -> Murder Dym Dxm DyM DxM * ── 3. Cluster-robust SEs ───────────────────────────────────────── * Stata's vce(cluster state) on -regress- applies exactly the * small-sample correction used in the Fitzgerald et al. (2026) * replication code: * * V_cluster = (N-1)/(N-k) * G/(G-1) * (X'X)^-1 . S . (X'X)^-1 * * with G = 48 clusters of 12 observations each. No extra code needed. * ── 4. Estimator A - First-difference OLS (no controls) ────────── * ESTIMAND. Throughout the five estimators below the parameter of * interest is alpha, the average partial effect of (first-differenced) * effective abortion rate on (first-differenced) state crime rate. * Identification rests on: * (1) conditional independence given the 284 partialled controls; and * (2) parallel trends in levels (state FEs absorbed by differencing, * year FEs by the partialling step). di _n(2) "========================================" di "STEP 4 - FIRST-DIFFERENCE OLS (no controls)" di "========================================" * Storage for downstream table. tempname FD_b FD_se mat `FD_b' = J(3, 1, .) mat `FD_se' = J(3, 1, .) local row = 0 foreach o in v p m { local ++row local lab = cond("`o'"=="v","Violent crime",cond("`o'"=="p","Property crime","Murder")) local Y = cond("`o'"=="v","Dyv", cond("`o'"=="p","Dyp","Dym")) local D = cond("`o'"=="v","Dxv", cond("`o'"=="p","Dxp","Dxm")) di _n " Outcome: `lab' (y=`Y', d=`D')" qui regress `Y' `D', noconstant vce(cluster state) mat `FD_b'[`row',1] = _b[`D'] mat `FD_se'[`row',1] = _se[`D'] di as text " alpha_hat = " as result %9.4f _b[`D'] /// as text " SE = " as result %9.4f _se[`D'] } * ── 5. Estimator B - OLS with all 284 controls ─────────────────── di _n(2) "========================================" di "STEP 5 - OLS WITH ALL 284 CONTROLS (kitchen-sink)" di "========================================" di "Feasible because p=284 < n=576: OLS technically inverts. But" di "many controls are near-collinear, so SEs balloon. This is what" di "motivates LASSO: keep the controls that matter, drop the rest." tempname OLS_b OLS_se OLS_n mat `OLS_b' = J(3, 1, .) mat `OLS_se' = J(3, 1, .) mat `OLS_n' = J(3, 1, .) local row = 0 foreach o in v p m { local ++row local lab = cond("`o'"=="v","Violent crime",cond("`o'"=="p","Property crime","Murder")) local Y = cond("`o'"=="v","DyV", cond("`o'"=="p","DyP","DyM")) local D = cond("`o'"=="v","DxV", cond("`o'"=="p","DxP","DxM")) di _n " Outcome: `lab' (y=`Y', d=`D')" qui regress `Y' `D' z`o'1-z`o'284, noconstant vce(cluster state) mat `OLS_b'[`row',1] = _b[`D'] mat `OLS_se'[`row',1] = _se[`D'] mat `OLS_n'[`row',1] = 284 di as text " alpha_hat = " as result %9.4f _b[`D'] /// as text " SE = " as result %9.4f _se[`D'] /// as text " (using 284 controls)" } * ── 6. Estimator C - Post-Structural LASSO (PSL) ───────────────── * PSL = one LASSO on (d, X) -> y with the treatment d FORCED IN * (penalty.factor=0 in R; pnotpen() in Stata's rlasso). The selected * controls are passed to plain OLS for the final alpha. * * DESIGN NOTE: The R companion uses cv.glmnet for PSL (3-fold CV). * Stata's cvlasso supports notpen(), but its CV path through 100+ * lambda values with notpen-partialling is dramatically slower than * cv.glmnet's coordinate descent. We use the Belloni-Chernozhukov- * Hansen RIGOROUS penalty via rlasso instead, with d pinned via * pnotpen(). This is fast, deterministic, and answers the same * conceptual question PSL was designed to ask: "one LASSO with the * treatment forced in, then post-OLS". The penalty rule differs * across the two implementations but the recipe is identical. * * WHY POST-OLS, not LASSO coefficients? LASSO shrinks the * coefficients of variables it keeps toward zero - that introduces * bias in alpha. Refitting with plain OLS on the selected variables * removes the shrinkage. Throughout this script LASSO is used for * SELECTION only; the final alpha always comes from OLS. di _n(2) "========================================" di "STEP 6 - POST-STRUCTURAL LASSO (PSL)" di "========================================" di "One rlasso on (d, X) -> y with d pinned (pnotpen)," di "then OLS on d + selected controls. Rigorous penalty (c=1.1, gamma=0.05)." tempname PSL_b PSL_se PSL_n mat `PSL_b' = J(3, 1, .) mat `PSL_se' = J(3, 1, .) mat `PSL_n' = J(3, 1, .) local row = 0 foreach o in v p m { local ++row local lab = cond("`o'"=="v","Violent crime",cond("`o'"=="p","Property crime","Murder")) local Y = cond("`o'"=="v","DyV", cond("`o'"=="p","DyP","DyM")) local D = cond("`o'"=="v","DxV", cond("`o'"=="p","DxP","DxM")) di _n " Outcome: `lab' (y=`Y', d=`D')" qui rlasso `Y' `D' z`o'1-z`o'284, nocons pnotpen(`D') c(1.1) gamma(0.05) local sel "`e(selected)'" * e(selected) includes pnotpen variable(s) (here: DxV). Strip d out * so we can pass clean "d + controls" to the post-OLS regression. local sel : list sel - D local nsel : word count `sel' qui regress `Y' `D' `sel', noconstant vce(cluster state) mat `PSL_b'[`row',1] = _b[`D'] mat `PSL_se'[`row',1] = _se[`D'] mat `PSL_n'[`row',1] = `nsel' di as text " alpha_hat = " as result %9.4f _b[`D'] /// as text " SE = " as result %9.4f _se[`D'] /// as text " | `nsel' controls selected" } * ── 7. Estimator D - Double LASSO, rigorous penalty (rlasso) ───── * Belloni-Chernozhukov-Hansen Double LASSO with the rigorous penalty: * * 1. rlasso y on X -> I_y = selected outcome-equation indices * 2. rlasso d on X -> I_d = selected treatment-equation indices * 3. OLS y on d + X[I_y ∪ I_d] with state-clustered SEs * * The rigorous penalty is data-driven (Belloni, Chen, Chernozhukov & * Hansen 2012) and chosen so that selection-error noise is dominated * by the signal. It is much more parsimonious than CV. * * Penalty constants c=1.1 and gamma=0.05 match the JAE (2026) * replication code (readdata_all_OLS.R lines 585, 653) and the * R companion's hdm::rlasso call. di _n(2) "========================================" di "STEP 7 - DOUBLE LASSO, RIGOROUS PENALTY (rlasso)" di "========================================" di "Two rlasso calls (y on X, d on X), union of selected, then post-OLS." di "'Rigorous' = lambda from Belloni et al. (2012) theory, not CV." tempname DLR_b DLR_se DLR_n DLR_Iy DLR_Id DLR_U mat `DLR_b' = J(3, 1, .) mat `DLR_se' = J(3, 1, .) mat `DLR_n' = J(3, 1, .) mat `DLR_Iy' = J(3, 1, .) mat `DLR_Id' = J(3, 1, .) mat `DLR_U' = J(3, 1, .) local row = 0 foreach o in v p m { local ++row local lab = cond("`o'"=="v","Violent crime",cond("`o'"=="p","Property crime","Murder")) local Y = cond("`o'"=="v","DyV", cond("`o'"=="p","DyP","DyM")) local D = cond("`o'"=="v","DxV", cond("`o'"=="p","DxP","DxM")) di _n " Outcome: `lab' (y=`Y', d=`D')" * Step 1: LASSO y on X. qui rlasso `Y' z`o'1-z`o'284, nocons c(1.1) gamma(0.05) local Iy "`e(selected)'" local nIy : word count `Iy' * Step 2: LASSO d on X. qui rlasso `D' z`o'1-z`o'284, nocons c(1.1) gamma(0.05) local Id "`e(selected)'" local nId : word count `Id' * Step 3: union of selected, then post-OLS. local U : list Iy | Id local nU : word count `U' if `nU' > 0 { qui regress `Y' `D' `U', noconstant vce(cluster state) } else { * Fall back to univariate first-difference fit (no controls survived). qui regress `Y' `D', noconstant vce(cluster state) } mat `DLR_b'[`row',1] = _b[`D'] mat `DLR_se'[`row',1] = _se[`D'] mat `DLR_n'[`row',1] = `nU' mat `DLR_Iy'[`row',1] = `nIy' mat `DLR_Id'[`row',1] = `nId' mat `DLR_U'[`row',1] = `nU' di as text " |I_y| = " as result %3.0f `nIy' /// as text " |I_d| = " as result %3.0f `nId' /// as text " |union| = " as result %3.0f `nU' di as text " alpha_hat = " as result %9.4f _b[`D'] /// as text " SE = " as result %9.4f _se[`D'] } * ── 8. Estimator E - Double LASSO, CV penalty (cvlasso) ────────── * Same three steps as section 7, but each LASSO is tuned by 3-fold * CV (matching Fitzgerald et al. 2026 footnote 2). Lambda is picked * to minimise out-of-sample MSE (lopt = "lambda at MSE minimum"). * * RUNTIME NOTE. Stata's cvlasso at p/n ≈ 0.5 is dramatically slower * than R's cv.glmnet: a single call with the default lcount(100) * takes 5+ minutes per outcome-equation, which makes the 6-call DL-CV * pipeline impractical. We use lcount(10) to get the CV grid down to * 10 lambda values, sacrificing precision in lambda selection in * exchange for finishing in under 60 seconds per call. The Stata-vs-R * drift in Section 14 captures the consequence: DL-CV's α and selected * set sizes diverge more than the other Tier-C estimators do. di _n(2) "========================================" di "STEP 8 - DOUBLE LASSO, CV PENALTY (cvlasso, 3-fold, lcount=10)" di "========================================" tempname DLC_b DLC_se DLC_n DLC_Iy DLC_Id DLC_U mat `DLC_b' = J(3, 1, .) mat `DLC_se' = J(3, 1, .) mat `DLC_n' = J(3, 1, .) mat `DLC_Iy' = J(3, 1, .) mat `DLC_Id' = J(3, 1, .) mat `DLC_U' = J(3, 1, .) local row = 0 foreach o in v p m { local ++row local lab = cond("`o'"=="v","Violent crime",cond("`o'"=="p","Property crime","Murder")) local Y = cond("`o'"=="v","DyV", cond("`o'"=="p","DyP","DyM")) local D = cond("`o'"=="v","DxV", cond("`o'"=="p","DxP","DxM")) di _n " Outcome: `lab' (y=`Y', d=`D')" qui cvlasso `Y' z`o'1-z`o'284, nfolds(3) seed(20260520) lopt lglmnet lcount(10) local Iy "`e(selected)'" local nIy : word count `Iy' qui cvlasso `D' z`o'1-z`o'284, nfolds(3) seed(20260520) lopt lglmnet lcount(10) local Id "`e(selected)'" local nId : word count `Id' local U : list Iy | Id local nU : word count `U' if `nU' > 0 { qui regress `Y' `D' `U', noconstant vce(cluster state) } else { qui regress `Y' `D', noconstant vce(cluster state) } mat `DLC_b'[`row',1] = _b[`D'] mat `DLC_se'[`row',1] = _se[`D'] mat `DLC_n'[`row',1] = `nU' mat `DLC_Iy'[`row',1] = `nIy' mat `DLC_Id'[`row',1] = `nId' mat `DLC_U'[`row',1] = `nU' di as text " |I_y| = " as result %3.0f `nIy' /// as text " |I_d| = " as result %3.0f `nId' /// as text " |union| = " as result %3.0f `nU' di as text " alpha_hat = " as result %9.4f _b[`D'] /// as text " SE = " as result %9.4f _se[`D'] } * ── 9. Build results_table2.csv (5 estimators x 3 outcomes) ────── di _n(2) "========================================" di "STEP 9 - REPLICATION OF PAPER TABLE 2" di "========================================" preserve clear set obs 15 gen str14 method = "" gen str16 outcome = "" gen double estimate = . gen double std_error = . gen long n_selected = . local outcomes "viol prop murd" local methods "FD OLS PSL DLR DLC" local mlabels `""First diff" "OLS (full)" "PSL" "DL (rigorous)" "DL (CV)""' * Fill the table by hand, row-by-row, looking up the right matrix per cell. local i = 0 forvalues oi = 1/3 { local o : word `oi' of `outcomes' if "`o'" == "viol" local olab "Violent crime" if "`o'" == "prop" local olab "Property crime" if "`o'" == "murd" local olab "Murder" forvalues mi = 1/5 { local ++i local m : word `mi' of `methods' local mlab : word `mi' of `mlabels' qui replace method = "`mlab'" in `i' qui replace outcome = "`olab'" in `i' if "`m'" == "FD" { qui replace estimate = `FD_b'[`oi',1] in `i' qui replace std_error = `FD_se'[`oi',1] in `i' qui replace n_selected = 0 in `i' } else if "`m'" == "OLS" { qui replace estimate = `OLS_b'[`oi',1] in `i' qui replace std_error = `OLS_se'[`oi',1] in `i' qui replace n_selected = `OLS_n'[`oi',1] in `i' } else if "`m'" == "PSL" { qui replace estimate = `PSL_b'[`oi',1] in `i' qui replace std_error = `PSL_se'[`oi',1] in `i' qui replace n_selected = `PSL_n'[`oi',1] in `i' } else if "`m'" == "DLR" { qui replace estimate = `DLR_b'[`oi',1] in `i' qui replace std_error = `DLR_se'[`oi',1] in `i' qui replace n_selected = `DLR_n'[`oi',1] in `i' } else if "`m'" == "DLC" { qui replace estimate = `DLC_b'[`oi',1] in `i' qui replace std_error = `DLC_se'[`oi',1] in `i' qui replace n_selected = `DLC_n'[`oi',1] in `i' } } } gen double ci_lo = estimate - 1.96 * std_error gen double ci_hi = estimate + 1.96 * std_error order method outcome estimate std_error n_selected ci_lo ci_hi list, sepby(outcome) abbreviate(20) export delimited "results_table2.csv", replace di _n "Wrote results_table2.csv" restore * ── 9b. Selection diagnostic CSV ───────────────────────────────── preserve clear set obs 6 gen str16 outcome = "" gen str14 method = "" gen long n_Iy = . gen long n_Id = . gen long n_union = . local i = 0 forvalues oi = 1/3 { local o : word `oi' of viol prop murd if "`o'" == "viol" local olab "Violent crime" if "`o'" == "prop" local olab "Property crime" if "`o'" == "murd" local olab "Murder" local ++i qui replace outcome = "`olab'" in `i' qui replace method = "DL (rigorous)" in `i' qui replace n_Iy = `DLR_Iy'[`oi',1] in `i' qui replace n_Id = `DLR_Id'[`oi',1] in `i' qui replace n_union = `DLR_U'[`oi',1] in `i' local ++i qui replace outcome = "`olab'" in `i' qui replace method = "DL (CV)" in `i' qui replace n_Iy = `DLC_Iy'[`oi',1] in `i' qui replace n_Id = `DLC_Id'[`oi',1] in `i' qui replace n_union = `DLC_U'[`oi',1] in `i' } list, sepby(outcome) abbreviate(20) export delimited "selection_diagnostic.csv", replace di "Wrote selection_diagnostic.csv" restore * ── 10. Figures ────────────────────────────────────────────────── * * Figure construction is delegated to the companion file figures.do, * which reads the two CSVs we just wrote (results_table2.csv and * selection_diagnostic.csv) and rebuilds all three PNGs in seconds. * * Keeping the plotting logic in one file gives us a single source of * truth: a styling tweak only needs to be made in figures.do (no need * to re-run the full ~10-minute LASSO pipeline). * * The three outputs are: * stata_double_lasso_estimates.png (forest plot, 3 panels) * stata_double_lasso_selection.png (selection bars, 4 metrics x 2 methods) * stata_double_lasso_methods_compare.png (rigorous vs CV with point labels) di _n(2) "========================================" di "STEP 10 - FIGURES (delegated to figures.do)" di "========================================" do figures.do * ── 11. Summary ────────────────────────────────────────────────── di _n(2) "========================================" di "STEP 11 - SUMMARY" di "========================================" di "" di "Generated files:" di " results_table2.csv (5 estimators x 3 outcomes)" di " selection_diagnostic.csv (|I_y|, |I_d|, union)" di " stata_double_lasso_estimates.png (forest plot)" di " stata_double_lasso_selection.png (selection bars)" di " stata_double_lasso_methods_compare.png (rigorous vs CV)" di "" di "Key practical insight from Fitzgerald et al. (2026):" di "DL helps most when the TREATMENT is highly predictable from the" di "controls but the OUTCOME is not. That is the case here: the" di "effective abortion rate is well explained by lagged demographics" di "and within-state trends, while crime is much noisier." * ── 12. pdslasso one-line demo (three-estimator output) ────────── * The explicit DL-rigorous recipe in Step 7 corresponds to the PDS * (post-double-selection) column that pdslasso reports. But pdslasso * also computes two additional valid estimators — the lasso- * orthogonalized and post-lasso-orthogonalized versions — based on * the Belloni-Chernozhukov-Hansen-Chen orthogonalization framework. * All three target the same causal alpha; they differ in how the * high-dimensional controls are residualised out. The post's Section * 8 walks through the three side-by-side; the call below produces * the actual three-panel output on the violent-crime equation. di _n(2) "========================================" di "STEP 12 - pdslasso three-estimator demo (violent crime)" di "========================================" pdslasso DyV DxV (zv1-zv284), cluster(state) loptions(c(1.1) gamma(0.05)) di _n "=== Script completed successfully ===" log close