------------------------------------------------------------------------------- name: log: /Users/carlosmendez/Documents/Github/starter-academic-v501/content > /post/stata_bma_dsl/analysis.log log type: text opened on: 31 Mar 2026, 15:38:50 . . display _newline "=============================================" ============================================= . display " TUTORIAL: BMA and DSL for the EKC" TUTORIAL: BMA and DSL for the EKC . display " Data: Synthetic panel (80 countries, 1995-2014)" Data: Synthetic panel (80 countries, 1995-2014) . display " Started: $S_DATE $S_TIME" Started: 31 Mar 2026 15:38:50 . display "=============================================" ============================================= . . * Load synthetic data from GitHub for reproducibility . import delimited "https://github.com/cmg777/starter-academic-v501/raw/master/ > content/post/stata_bma_dsl/synthetic_ekc_panel.csv", clear (encoding automatically selected: ISO-8859-1) (18 vars, 1,600 obs) . . xtset country_id year, yearly Panel variable: country_id (strongly balanced) Time variable: year, 1995 to 2014 Delta: 1 year . . * Label variables for readable output . label variable country_id "Country ID" . label variable year "Year" . label variable ln_co2 "CO2 per capita (log)" . label variable ln_gdp "GDP per capita (log)" . label variable ln_gdp_sq "GDP per capita squared (log)" . label variable ln_gdp_cb "GDP per capita cubed (log)" . label variable fossil_fuel "Fossil fuel share (%)" . label variable renewable "Renewable energy (%)" . label variable urban "Urban population (%)" . label variable globalization "Globalization index" . label variable pop_density "Population density" . label variable democracy "Democracy score" . label variable corruption "Corruption index" . label variable industry "Industry VA (% GDP)" . label variable services "Services VA (% GDP)" . label variable trade "Trade openness (% GDP)" . label variable fdi "FDI inflows (% GDP)" . label variable credit "Domestic credit (% GDP)" . . display _newline "=== Dataset summary ===" === Dataset summary === . xtdescribe country_id: 1, 2, ..., 80 n = 80 year: 1995, 1996, ..., 2014 T = 20 Delta(year) = 1 year Span(year) = 20 periods (country_id*year uniquely identifies each observation) Distribution of T_i: min 5% 25% 50% 75% 95% max 20 20 20 20 20 20 20 Freq. Percent Cum. | Pattern ---------------------------+---------------------- 80 100.00 100.00 | 11111111111111111111 ---------------------------+---------------------- 80 100.00 | XXXXXXXXXXXXXXXXXXXX . summarize $outcome $gdp_vars $controls Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- ln_co2 | 1,600 -19.0385 .7863276 -21.03685 -16.8315 ln_gdp | 1,600 9.583865 1.329675 6.974263 11.9704 ln_gdp_sq | 1,600 93.61741 25.55106 48.64035 143.2904 ln_gdp_cb | 1,600 931.1045 373.829 339.2306 1715.243 fossil_fuel | 1,600 54.7724 19.14168 6.36807 95 -------------+--------------------------------------------------------- renewable | 1,600 29.54128 11.96568 1 64.2207 urban | 1,600 53.67416 14.778 15.95174 91.63234 globalizat~n | 1,600 57.64978 12.71537 26.75758 95 pop_density | 1,600 121.3443 210.2646 1 1571.771 democracy | 1,600 2.333462 4.179503 -6.12244 10 -------------+--------------------------------------------------------- corruption | 1,600 52.35226 28.52792 0 100 industry | 1,600 24.64331 6.180478 5.843938 45.32926 services | 1,600 43.55976 9.366089 17.82623 64.07455 trade | 1,600 67.43547 19.36148 10.04306 128.0595 fdi | 1,600 2.982371 4.373857 -11.50437 16.19903 -------------+--------------------------------------------------------- credit | 1,600 53.44019 18.20204 11.32991 123.2399 . . . *============================================================================ > =* . * SECTION 2: EXPLORE THE DATA -- FIGURE 1 . *============================================================================ > =* . . twoway (scatter $outcome ln_gdp, /// > msize(vsmall) mcolor("106 155 204"%40) msymbol(circle)), /// > ytitle("Log CO2 per capita") /// > xtitle("Log GDP per capita") /// > title("Synthetic Data: CO2 vs. Income", size(medium)) /// > subtitle("80 countries, 1995-2014 (N = 1,600)", size(small)) /// > note("Each dot is a country-year observation." /// > "The nonlinear pattern motivates testing for a cubic EKC shape.") // > / > scheme(s2color) /// > name(fig1_scatter, replace) . . graph export "stata_bma_dsl_fig1_scatter.png", replace width(2400) file stata_bma_dsl_fig1_scatter.png written in PNG format . display _newline "Saved: stata_bma_dsl_fig1_scatter.png" Saved: stata_bma_dsl_fig1_scatter.png . . . *============================================================================ > =* . * SECTION 3: BASELINE FIXED EFFECTS . *============================================================================ > =* . . display _newline "=============================================" ============================================= . display " BASELINE: Standard Fixed Effects" BASELINE: Standard Fixed Effects . display "=============================================" ============================================= . . *---------------------------------------------* . * 3a. Sparse specification (no controls) * . *---------------------------------------------* . reghdfe $outcome $gdp_vars, absorb(country_id year) vce(cluster country_id) (MWFE estimator converged in 2 iterations) HDFE Linear regression Number of obs = 1,600 Absorbing 2 HDFE groups F( 3, 79) = 16.62 Statistics robust to heteroskedasticity Prob > F = 0.0000 R-squared = 0.9620 Adj R-squared = 0.9594 Within R-sq. = 0.0354 Number of clusters (country_id) = 80 Root MSE = 0.1584 (Std. err. adjusted for 80 clusters in country_id) ------------------------------------------------------------------------------ | Robust ln_co2 | Coefficient std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- ln_gdp | -7.498046 1.623988 -4.62 0.000 -10.73051 -4.26558 ln_gdp_sq | .848967 .1704533 4.98 0.000 .5096881 1.188246 ln_gdp_cb | -.0314993 .005931 -5.31 0.000 -.0433047 -.019694 _cons | 2.672866 5.155102 0.52 0.606 -7.588108 12.93384 ------------------------------------------------------------------------------ Absorbed degrees of freedom: -----------------------------------------------------+ Absorbed FE | Categories - Redundant = Num. Coefs | -------------+---------------------------------------| country_id | 80 80 0 *| year | 20 1 19 | -----------------------------------------------------+ * = FE nested within cluster; treated as redundant for DoF computation . estimates store fe_sparse . . local b1_sparse = _b[ln_gdp] . local b2_sparse = _b[ln_gdp_sq] . local b3_sparse = _b[ln_gdp_cb] . . display _newline "=== Sparse FE Coefficients ===" === Sparse FE Coefficients === . display " b1 (GDP): " %9.4f `b1_sparse' b1 (GDP): -7.4980 . display " b2 (GDP^2): " %9.4f `b2_sparse' b2 (GDP^2): 0.8490 . display " b3 (GDP^3): " %9.4f `b3_sparse' b3 (GDP^3): -0.0315 . . *---------------------------------------------* . * 3b. Kitchen-sink specification * . *---------------------------------------------* . reghdfe $outcome $gdp_vars $controls, absorb(country_id year) vce(cluster cou > ntry_id) (MWFE estimator converged in 2 iterations) HDFE Linear regression Number of obs = 1,600 Absorbing 2 HDFE groups F( 15, 79) = 15.77 Statistics robust to heteroskedasticity Prob > F = 0.0000 R-squared = 0.9655 Adj R-squared = 0.9629 Within R-sq. = 0.1249 Number of clusters (country_id) = 80 Root MSE = 0.1515 (Std. err. adjusted for 80 clusters in country_id) ------------------------------------------------------------------------------ | Robust ln_co2 | Coefficient std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- ln_gdp | -7.130693 1.562581 -4.56 0.000 -10.24093 -4.020453 ln_gdp_sq | .8059928 .1647973 4.89 0.000 .477972 1.134014 ln_gdp_cb | -.0298133 .0057365 -5.20 0.000 -.0412314 -.0183951 fossil_fuel | .0138444 .0014853 9.32 0.000 .010888 .0168008 renewable | -.006795 .0019322 -3.52 0.001 -.0106409 -.0029491 urban | .0057534 .0021432 2.68 0.009 .0014875 .0100192 globalizat~n | .0015186 .0012832 1.18 0.240 -.0010357 .0040728 pop_density | .0000794 .0002303 0.34 0.731 -.000379 .0005378 democracy | -.0002971 .007735 -0.04 0.969 -.0156933 .0150991 corruption | .0009812 .0008415 1.17 0.247 -.0006936 .0026561 industry | .0086336 .0017848 4.84 0.000 .0050811 .0121861 services | -.0005642 .0017205 -0.33 0.744 -.0039889 .0028604 trade | -.0002458 .0007695 -0.32 0.750 -.0017774 .0012858 fdi | -.0017599 .0019509 -0.90 0.370 -.005643 .0021232 credit | -.00139 .0007516 -1.85 0.068 -.002886 .0001061 _cons | .499146 4.926364 0.10 0.920 -9.306537 10.30483 ------------------------------------------------------------------------------ Absorbed degrees of freedom: -----------------------------------------------------+ Absorbed FE | Categories - Redundant = Num. Coefs | -------------+---------------------------------------| country_id | 80 80 0 *| year | 20 1 19 | -----------------------------------------------------+ * = FE nested within cluster; treated as redundant for DoF computation . estimates store fe_kitchen . . local b1_kitchen = _b[ln_gdp] . local b2_kitchen = _b[ln_gdp_sq] . local b3_kitchen = _b[ln_gdp_cb] . . display _newline "=== Kitchen-Sink FE Coefficients ===" === Kitchen-Sink FE Coefficients === . display " b1 (GDP): " %9.4f `b1_kitchen' b1 (GDP): -7.1307 . display " b2 (GDP^2): " %9.4f `b2_kitchen' b2 (GDP^2): 0.8060 . display " b3 (GDP^3): " %9.4f `b3_kitchen' b3 (GDP^3): -0.0298 . . *---------------------------------------------* . * 3c. Compare specifications * . *---------------------------------------------* . display _newline "=== Coefficient Comparison ===" === Coefficient Comparison === . display _newline " Sparse FE Kitchen-Sink FE" Sparse FE Kitchen-Sink FE . display " b1 (GDP): " %9.4f `b1_sparse' " " %9.4f `b1_kitchen' b1 (GDP): -7.4980 -7.1307 . display " b2 (GDP^2): " %9.4f `b2_sparse' " " %9.4f `b2_kitchen' b2 (GDP^2): 0.8490 0.8060 . display " b3 (GDP^3): " %9.4f `b3_sparse' " " %9.4f `b3_kitchen' b3 (GDP^3): -0.0315 -0.0298 . . * Turning points . local disc_sparse = (`b2_sparse')^2 - 3 * `b1_sparse' * `b3_sparse' . local disc_kitchen = (`b2_kitchen')^2 - 3 * `b1_kitchen' * `b3_kitchen' . . display _newline "=== Turning Points ===" === Turning Points === . if `disc_sparse' > 0 { . local min_sparse = exp((-`b2_sparse' + sqrt(`disc_sparse')) / (3 * `b3_sp > arse')) . local max_sparse = exp((-`b2_sparse' - sqrt(`disc_sparse')) / (3 * `b3_sp > arse')) . display " Sparse: Min $" %8.0fc `min_sparse' " Max $" %8.0fc `max_spar > se' Sparse: Min $ 2,478 Max $ 25,656 . } . else { . display " Sparse: No real turning points" . local min_sparse = . . local max_sparse = . . } . if `disc_kitchen' > 0 { . local min_kitchen = exp((-`b2_kitchen' + sqrt(`disc_kitchen')) / (3 * `b3 > _kitchen')) . local max_kitchen = exp((-`b2_kitchen' - sqrt(`disc_kitchen')) / (3 * `b3 > _kitchen')) . display " Kitchen: Min $" %8.0fc `min_kitchen' " Max $" %8.0fc `max_kit > chen' Kitchen: Min $ 2,426 Max $ 27,694 . } . else { . display " Kitchen: No real turning points" . local min_kitchen = . . local max_kitchen = . . } . . *---------------------------------------------* . * 3d. Coefficient instability chart -- Fig 2 * . *---------------------------------------------* . preserve . clear . set obs 6 Number of observations (_N) was 0, now 6. . . gen spec = "" (6 missing values generated) . gen coef_name = "" (6 missing values generated) . gen value = . (6 missing values generated) . gen order = _n . . replace spec = "Sparse FE" in 1 variable spec was str1 now str9 (1 real change made) . replace coef_name = "b1 (GDP)" in 1 variable coef_name was str1 now str8 (1 real change made) . replace value = `b1_sparse' in 1 (1 real change made) . . replace spec = "Sparse FE" in 2 (1 real change made) . replace coef_name = "b2 (GDP sq)" in 2 variable coef_name was str8 now str11 (1 real change made) . replace value = `b2_sparse' in 2 (1 real change made) . . replace spec = "Sparse FE" in 3 (1 real change made) . replace coef_name = "b3 (GDP cb)" in 3 (1 real change made) . replace value = `b3_sparse' in 3 (1 real change made) . . replace spec = "Kitchen-Sink FE" in 4 variable spec was str9 now str15 (1 real change made) . replace coef_name = "b1 (GDP)" in 4 (1 real change made) . replace value = `b1_kitchen' in 4 (1 real change made) . . replace spec = "Kitchen-Sink FE" in 5 (1 real change made) . replace coef_name = "b2 (GDP sq)" in 5 (1 real change made) . replace value = `b2_kitchen' in 5 (1 real change made) . . replace spec = "Kitchen-Sink FE" in 6 (1 real change made) . replace coef_name = "b3 (GDP cb)" in 6 (1 real change made) . replace value = `b3_kitchen' in 6 (1 real change made) . . graph twoway /// > (bar value order if spec == "Sparse FE", /// > barwidth(0.35) color("106 155 204")) /// > (bar value order if spec == "Kitchen-Sink FE", /// > barwidth(0.35) color("217 119 87")), /// > xlabel(1 `""b1" "(GDP)""' 2 `""b2" "(GDP{sup:2})""' 3 `""b3" "(GDP{sup:3} > )""' /// > 4 `""b1" "(GDP)""' 5 `""b2" "(GDP{sup:2})""' 6 `""b3" "(GDP{sup:3} > )""', /// > angle(0) labsize(small)) /// > xline(3.5, lcolor(gs10) lpattern(dash)) /// > ylabel(, format(%5.2f)) /// > ytitle("Coefficient value") /// > title("Coefficient Instability Across Specifications", size(medium)) /// > subtitle("Same GDP terms, different control sets", size(small)) /// > legend(order(1 "Sparse FE (no controls)" 2 "Kitchen-Sink FE (all 12 contr > ols)") /// > rows(1) position(6) size(small)) /// > note("Adding controls shifts GDP coefficients --- which specification is > correct?") /// > scheme(s2color) /// > name(fig2_instab, replace) . . graph export "stata_bma_dsl_fig2_instability.png", replace width(2400) file stata_bma_dsl_fig2_instability.png written in PNG format . display _newline "Saved: stata_bma_dsl_fig2_instability.png" Saved: stata_bma_dsl_fig2_instability.png . restore . . . *============================================================================ > =* . * SECTION 4: BAYESIAN MODEL AVERAGING (BMA) . *============================================================================ > =* . . display _newline "=============================================" ============================================= . display " BMA ESTIMATION" BMA ESTIMATION . display " Starting: $S_DATE $S_TIME" Starting: 31 Mar 2026 15:38:53 . display "=============================================" ============================================= . . bmaregress $outcome $gdp_vars $controls /// > ($fe, always), /// > mprior(uniform) groupfv gprior(uip) /// > mcmcsize(50000) dots(5000, every(10000)) /// > rseed(9988) inputorder pipcutoff(0.8) /// > saving("_bma_temp.dta", replace) Burn-in (2500): done Simulation (50000): .10000.20000.30000.40000.50000 done Computing model probabilities ... Bayesian model averaging No. of obs = 1,600 Linear regression No. of predictors = 113 MC3 sampling Groups = 15 Always = 98 No. of models = 163 For CPMP >= .9 = 14 Priors: Mean model size = 104.578 Models: Uniform Burn-in = 2,500 Cons.: Noninformative MCMC sample size = 50,000 Coef.: Zellner's g Acceptance rate = 0.0904 g: Unit-information, g = 1,600 Shrinkage, g/(1+g) = 0.9994 sigma2: Noninformative Mean sigma2 = 0.022 Sampling correlation = 0.9997 ------------------------------------------------------------------------------ ln_co2 | Mean Std. dev. Group PIP -------------+---------------------------------------------------------------- ln_gdp | -7.13901 1.811093 1 .99401 ln_gdp_sq | .8078437 .1892418 2 .99991 ln_gdp_cb | -.0299182 .0065105 3 .99976 fossil_fuel | .0138139 .001283 4 1 renewable | -.0068332 .0023506 5 .95945 industry | .0085503 .0019766 11 .99867 -------------+---------------------------------------------------------------- Always | country_id | 2 | -.1179799 .1145024 0 1 3 | -.034187 .112363 0 1 4 | -.6996246 .201126 0 1 5 | .829754 .0962586 0 1 6 | -.0613427 .1977392 0 1 7 | -.106526 .1839812 0 1 8 | -.7705366 .1379249 0 1 9 | -.2429933 .0680603 0 1 10 | .5178117 .1919635 0 1 11 | -.7845254 .0788075 0 1 12 | .1495013 .1486364 0 1 13 | -.3565051 .1804139 0 1 14 | -.2302201 .1129516 0 1 15 | -1.39671 .0789442 0 1 16 | -.1258029 .119203 0 1 17 | -.9472684 .080517 0 1 18 | -.5984079 .1323325 0 1 19 | .0637678 .1512518 0 1 20 | .047098 .0749263 0 1 21 | -.7812588 .0754839 0 1 22 | -1.067352 .0694149 0 1 23 | -.4488343 .1185606 0 1 24 | -.1361108 .1251686 0 1 25 | -.7238642 .0923753 0 1 26 | .2767268 .1140441 0 1 27 | -.4335896 .1797908 0 1 28 | .0344115 .2160191 0 1 29 | -.8559822 .0736341 0 1 30 | -.6669463 .0588199 0 1 31 | -1.063821 .0735287 0 1 32 | -.0100271 .0715373 0 1 33 | -.0187314 .1794103 0 1 34 | -1.441324 .081423 0 1 35 | -.9283798 .0670747 0 1 36 | -.7646781 .0827938 0 1 37 | 1.124507 .178476 0 1 38 | -.7379022 .0603345 0 1 39 | -.5532597 .0691558 0 1 40 | .1730358 .1658478 0 1 41 | -.0445621 .1416842 0 1 42 | -.0645891 .0871661 0 1 43 | .0542686 .1401205 0 1 44 | .2162265 .1089101 0 1 45 | .2122891 .0646419 0 1 46 | .0155558 .0630781 0 1 47 | -.6256494 .0867368 0 1 48 | -.7090819 .1312863 0 1 49 | .3528996 .0727821 0 1 50 | -.6129385 .0933364 0 1 51 | -.7249113 .1131649 0 1 52 | .3607882 .1539638 0 1 53 | -.5461417 .0778006 0 1 54 | -1.386536 .071729 0 1 55 | -.8530744 .1779728 0 1 56 | -.0452429 .2174362 0 1 57 | -.6174655 .1009132 0 1 58 | -.5480472 .0728281 0 1 59 | .1442367 .0903137 0 1 60 | -.0574758 .1607972 0 1 61 | .7003917 .1163897 0 1 62 | .4690104 .187726 0 1 63 | .213479 .0663604 0 1 64 | -.5242577 .1286337 0 1 65 | -1.008797 .0644648 0 1 66 | -.2275966 .065082 0 1 67 | .6036724 .2213786 0 1 68 | -1.20648 .0676923 0 1 69 | -.2588551 .1749428 0 1 70 | -.0976466 .1188324 0 1 71 | -.2397095 .1954438 0 1 72 | .4040474 .1056883 0 1 73 | .0335636 .0835823 0 1 74 | -.3093181 .1042782 0 1 75 | -.1423568 .065284 0 1 76 | -.9109642 .0855667 0 1 77 | -.3019021 .2164476 0 1 78 | -.7025789 .1919114 0 1 79 | .236688 .1904104 0 1 80 | -.1709559 .1039995 0 1 | year | 1996 | -.0023553 .0234941 0 1 1997 | -.0016326 .0236192 0 1 1998 | .0057287 .0241149 0 1 1999 | .0109971 .0243942 0 1 2000 | -.009301 .0251807 0 1 2001 | -.0078331 .0256468 0 1 2002 | -.0283858 .0267905 0 1 2003 | -.0164654 .0273342 0 1 2004 | -.0461652 .0282498 0 1 2005 | -.0467576 .0293353 0 1 2006 | -.0446844 .0300448 0 1 2007 | -.0507042 .0316104 0 1 2008 | -.0759235 .0329535 0 1 2009 | -.0892308 .0338976 0 1 2010 | -.0889439 .0349373 0 1 2011 | -.0694394 .0364631 0 1 2012 | -.0976314 .0383168 0 1 2013 | -.1221188 .0389121 0 1 2014 | -.1000245 .0401478 0 1 | _cons | 1.093558 5.730563 0 1 ------------------------------------------------------------------------------ Note: Coefficient posterior means and std. dev. estimated from 163 models. Note: 9 predictors with PIP less than .8 not shown. Note: Acceptance rate is low. file _bma_temp.dta saved. . . display _newline "BMA completed: $S_DATE $S_TIME" BMA completed: 31 Mar 2026 15:39:05 . estimates store bma_uip . . *---------------------------------------------* . * 4a. Extract BMA coefficients * . *---------------------------------------------* . matrix bma_coefs = e(b_bma) . . local b1_bma = bma_coefs[1,1] . local b2_bma = bma_coefs[1,2] . local b3_bma = bma_coefs[1,3] . . display _newline "=== BMA Coefficient Signs ===" === BMA Coefficient Signs === . display " b1 (GDP): " %9.4f `b1_bma' _col(35) cond(`b1_bma' < 0, "NEGATIVE > ", "positive") b1 (GDP): -7.1390 NEGATIVE . display " b2 (GDP^2): " %9.4f `b2_bma' _col(35) cond(`b2_bma' > 0, "POSITIVE > ", "negative") b2 (GDP^2): 0.8078 POSITIVE . display " b3 (GDP^3): " %9.4f `b3_bma' _col(35) cond(`b3_bma' < 0, "NEGATIVE > ", "positive") b3 (GDP^3): -0.0299 NEGATIVE . . matrix bma_se = e(se_bma) . local se1_bma = bma_se[1,1] . local se2_bma = bma_se[1,2] . local se3_bma = bma_se[1,3] . . *---------------------------------------------* . * 4b. BMA turning points * . *---------------------------------------------* . local disc_bma = (`b2_bma')^2 - 3 * `b1_bma' * `b3_bma' . . display _newline "=== BMA Turning Points ===" === BMA Turning Points === . display " Discriminant: " %8.4f `disc_bma' /// > cond(`disc_bma' > 0, " (turning points exist)", " (no turning points)") Discriminant: 0.0119 (turning points exist) . . local min_bma = exp((-`b2_bma' + sqrt(`disc_bma')) / (3 * `b3_bma')) . local max_bma = exp((-`b2_bma' - sqrt(`disc_bma')) / (3 * `b3_bma')) . . display " Minimum: $" %8.0fc `min_bma' Minimum: $ 2,411 . display " Maximum: $" %8.0fc `max_bma' Maximum: $ 27,269 . . *---------------------------------------------* . * 4c. PIP bar chart -- Figure 3 * . *---------------------------------------------* . * Build a clean PIP chart with readable labels and color-coding . . matrix pip_mat = e(pip) . local nvars_all = colsof(pip_mat) . . preserve . clear . set obs `nvars_all' Number of observations (_N) was 0, now 116. . . gen varname = "" (116 missing values generated) . gen pip = . (116 missing values generated) . . local varnames : colnames pip_mat . forvalues i = 1/`nvars_all' { 2. local vname : word `i' of `varnames' 3. replace varname = "`vname'" in `i' 4. replace pip = pip_mat[1,`i'] in `i' 5. } variable varname was str1 now str6 (1 real change made) (1 real change made) variable varname was str6 now str9 (1 real change made) (1 real change made) (1 real change made) (1 real change made) variable varname was str9 now str11 (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) variable varname was str11 now str13 (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) . . * Keep ONLY the 15 candidate variables (drop FE dummies and constant) . keep if inlist(varname, "ln_gdp", "ln_gdp_sq", "ln_gdp_cb", /// > "fossil_fuel", "renewable", "urban", "globalization") | /// > inlist(varname, "pop_density", "democracy", "corruption", /// > "industry", "services", "trade", "fdi", "credit") (101 observations deleted) . local nvars = _N . . * Create readable labels . gen label = "" (15 missing values generated) . replace label = "GDP per capita (log)" if varname == "ln_gdp" variable label was str1 now str20 (1 real change made) . replace label = "GDP squared (log)" if varname == "ln_gdp_sq" (1 real change made) . replace label = "GDP cubed (log)" if varname == "ln_gdp_cb" (1 real change made) . replace label = "Fossil fuel share (%)" if varname == "fossil_fuel" variable label was str20 now str21 (1 real change made) . replace label = "Renewable energy (%)" if varname == "renewable" (1 real change made) . replace label = "Urban population (%)" if varname == "urban" (1 real change made) . replace label = "Globalization index" if varname == "globalization" (1 real change made) . replace label = "Population density" if varname == "pop_density" (1 real change made) . replace label = "Democracy score" if varname == "democracy" (1 real change made) . replace label = "Corruption index" if varname == "corruption" (1 real change made) . replace label = "Industry VA (% GDP)" if varname == "industry" (1 real change made) . replace label = "Services VA (% GDP)" if varname == "services" (1 real change made) . replace label = "Trade openness (% GDP)" if varname == "trade" variable label was str21 now str22 (1 real change made) . replace label = "FDI inflows (% GDP)" if varname == "fdi" (1 real change made) . replace label = "Domestic credit (% GDP)" if varname == "credit" variable label was str22 now str23 (1 real change made) . . * Mark true vs noise predictors . gen is_true = inlist(varname, "fossil_fuel", "renewable", "urban", /// > "democracy", "industry", "ln_gdp", "ln_gdp_sq", "ln_gdp_cb") . . * Sort by PIP (highest first) and create plotting order . gsort -pip . gen order = _n . . * Apply readable labels to the order variable . labmask order, values(label) . . * Display PIP values for the log . list varname label pip is_true, sep(0) +--------------------------------------------------------------+ | varname label pip is_true | |--------------------------------------------------------------| 1. | fossil_fuel Fossil fuel share (%) 1 1 | 2. | ln_gdp_sq GDP squared (log) .9999081 1 | 3. | ln_gdp_cb GDP cubed (log) .9997622 1 | 4. | industry Industry VA (% GDP) .9986705 1 | 5. | ln_gdp GDP per capita (log) .9940143 1 | 6. | renewable Renewable energy (%) .9594516 1 | 7. | urban Urban population (%) .2675386 1 | 8. | credit Domestic credit (% GDP) .1169078 0 | 9. | corruption Corruption index .0581883 0 | 10. | globalization Globalization index .0433192 0 | 11. | fdi FDI inflows (% GDP) .0410655 0 | 12. | services Services VA (% GDP) .0254143 0 | 13. | pop_density Population density .0253422 0 | 14. | trade Trade openness (% GDP) .0246842 0 | 15. | democracy Democracy score .0234042 1 | +--------------------------------------------------------------+ . . * Create the bar chart with color coding . graph twoway /// > (bar pip order if is_true == 1, horizontal barwidth(0.6) /// > color("106 155 204")) /// > (bar pip order if is_true == 0, horizontal barwidth(0.6) /// > color(gs11)), /// > xline(0.8, lcolor("217 119 87") lpattern(dash) lwidth(medium)) /// > ylabel(1(1)`nvars', valuelabel angle(0) labsize(small) nogrid) /// > xlabel(0(0.2)1, format(%3.1f)) /// > ytitle("") /// > xtitle("Posterior Inclusion Probability (PIP)") /// > title("BMA: Which Variables Matter?", size(medsmall)) /// > subtitle("Dashed line = 0.8 robustness threshold", size(small)) /// > legend(order(1 "True predictor (in DGP)" 2 "Noise variable (not in DGP)") > /// > rows(1) position(6) size(vsmall)) /// > note("Variables sorted by PIP. Blue = true predictor, gray = noise." /// > "Only candidate variables shown (country and year FE excluded).", si > ze(vsmall)) /// > scheme(s2color) ysize(7) xsize(9) /// > name(fig3_pip, replace) . . graph export "stata_bma_dsl_fig3_pip.png", replace width(2400) file stata_bma_dsl_fig3_pip.png written in PNG format . display _newline "Saved: stata_bma_dsl_fig3_pip.png" Saved: stata_bma_dsl_fig3_pip.png . restore . . *---------------------------------------------* . * 4d. Coefficient densities -- Figure 4 * . *---------------------------------------------* . * bmagraph coefdensity with multiple vars shows only the last one. . * Generate individual plots for all 6 robust variables (PIP > 0.80) . * and combine in a 3x2 grid. Use consistent small fonts, legends off. . . * All axis text at vsmall, zero reference line, no legend . local popts /// > xtitle("Coefficient", size(vsmall)) /// > ytitle("Density", size(vsmall)) /// > ylabel(, labsize(vsmall) angle(0) axis(1)) /// > ylabel(, labsize(vsmall) angle(0) axis(2)) /// > ytitle("Probability", size(vsmall) axis(2)) /// > xlabel(, labsize(vsmall)) /// > xline(0, lcolor(gs10) lpattern(dash) lwidth(thin)) /// > legend(off) scheme(s2color) . . bmagraph coefdensity ln_gdp, /// > title("GDP per capita (log)", size(small)) `popts' name(dens_gdp, replace > ) . bmagraph coefdensity ln_gdp_sq, /// > title("GDP squared (log)", size(small)) `popts' name(dens_gdp_sq, replace > ) . bmagraph coefdensity ln_gdp_cb, /// > title("GDP cubed (log)", size(small)) `popts' name(dens_gdp_cb, replace) . bmagraph coefdensity fossil_fuel, /// > title("Fossil fuel share (%)", size(small)) `popts' name(dens_fossil, rep > lace) . bmagraph coefdensity renewable, /// > title("Renewable energy (%)", size(small)) `popts' name(dens_renewable, r > eplace) . bmagraph coefdensity industry, /// > title("Industry VA (% GDP)", size(small)) `popts' name(dens_industry, rep > lace) . . graph combine dens_gdp dens_gdp_sq dens_gdp_cb /// > dens_fossil dens_renewable dens_industry, /// > cols(3) rows(2) imargin(small) /// > title("BMA: Posterior Coefficient Densities", size(medsmall)) /// > subtitle("All 6 robust variables (PIP > 0.80)", size(small)) /// > note("Blue curve = posterior density conditional on inclusion. Red line = > probability of noninclusion (1 - PIP)." /// > "Dashed gray line = zero. A blue curve far from zero = strong eviden > ce the variable matters.", size(vsmall)) /// > scheme(s2color) xsize(12) ysize(7) /// > name(fig4_density, replace) . . graph export "stata_bma_dsl_fig4_coefdensity.png", replace width(2400) file stata_bma_dsl_fig4_coefdensity.png written in PNG format . display _newline "Saved: stata_bma_dsl_fig4_coefdensity.png" Saved: stata_bma_dsl_fig4_coefdensity.png . . *---------------------------------------------* . * 4e. Pooled BMA (without fixed effects) * . *---------------------------------------------* . * Run BMA without country/year FE to show the cost of omitting FE. . * Without FE, noise variables get spuriously high PIPs. . . display _newline "=== BMA WITHOUT FIXED EFFECTS (POOLED) ===" === BMA WITHOUT FIXED EFFECTS (POOLED) === . bmaregress $outcome $gdp_vars $controls, /// > mprior(uniform) gprior(uip) /// > mcmcsize(50000) rseed(9988) pipcutoff(0.5) burnin(5000) Burn-in ... Simulation ... Computing model probabilities ... Bayesian model averaging No. of obs = 1,600 Linear regression No. of predictors = 15 MC3 sampling Groups = 15 Always = 0 No. of models = 36 For CPMP >= .9 = 7 Priors: Mean model size = 11.978 Models: Uniform Burn-in = 5,000 Cons.: Noninformative MCMC sample size = 50,000 Coef.: Zellner's g Acceptance rate = 0.0796 g: Unit-information, g = 1,600 Shrinkage, g/(1+g) = 0.9994 sigma2: Noninformative Mean sigma2 = 0.200 Sampling correlation = 0.9999 ------------------------------------------------------------------------------ ln_co2 | Mean Std. dev. Group PIP -------------+---------------------------------------------------------------- ln_gdp | -21.25798 1.641666 1 1 ln_gdp_sq | 2.284718 .1748828 2 1 ln_gdp_cb | -.0813933 .0061308 3 1 fossil_fuel | .0188852 .0010554 4 1 renewable | -.0192088 .0013911 5 1 urban | .0103139 .0012072 6 1 pop_density | -.0004314 .0000567 8 1 industry | .0138363 .0023477 11 1 services | .0164632 .0016573 12 1 credit | .0041023 .0008401 15 .99998 trade | -.0020936 .0010842 13 .85998 democracy | .0078792 .0042982 9 .84146 -------------+---------------------------------------------------------------- Always | _cons | 44.49613 5.11481 0 1 ------------------------------------------------------------------------------ Note: Coefficient posterior means and std. dev. estimated from 36 models. Note: 3 predictors with PIP less than .5 not shown. Note: Acceptance rate is low. . * Note: bmaregress does not support estimates store . matrix bma_pooled_tbl = r(table) . local b1_bma_p = bma_pooled_tbl[1,1] . local b2_bma_p = bma_pooled_tbl[1,2] . local b3_bma_p = bma_pooled_tbl[1,3] . local sd1_bma_p = bma_pooled_tbl[2,1] . local sd2_bma_p = bma_pooled_tbl[2,2] . local sd3_bma_p = bma_pooled_tbl[2,3] . . display _newline "Pooled BMA coefficients:" Pooled BMA coefficients: . display " b1 = " %9.4f `b1_bma_p' b1 = . . display " b2 = " %9.4f `b2_bma_p' b2 = . . display " b3 = " %9.4f `b3_bma_p' b3 = . . . * Turning points . local disc_bma_p = (`b2_bma_p')^2 - 3 * `b1_bma_p' * `b3_bma_p' . if `disc_bma_p' > 0 { . local min_bma_p = exp((-`b2_bma_p' + sqrt(`disc_bma_p')) / (3 * `b3_bma_p > ')) . local max_bma_p = exp((-`b2_bma_p' - sqrt(`disc_bma_p')) / (3 * `b3_bma_p > ')) . display " Minimum: $" %8.0fc `min_bma_p' Minimum: $ . . display " Maximum: $" %8.0fc `max_bma_p' Maximum: $ . . } . . . *============================================================================ > =* . * SECTION 5: DOUBLE-SELECTION LASSO (DSL) . *============================================================================ > =* . . display _newline "=============================================" ============================================= . display " DSL ESTIMATION" DSL ESTIMATION . display " Starting: $S_DATE $S_TIME" Starting: 31 Mar 2026 15:39:10 . display "=============================================" ============================================= . . dsregress $outcome $gdp_vars, /// > controls(($fe) $controls) /// > vce(cluster country_id) Estimating lasso for ln_co2 using plugin Estimating lasso for ln_gdp using plugin Estimating lasso for ln_gdp_sq using plugin Estimating lasso for ln_gdp_cb using plugin Double-selection linear model Number of obs = 1,600 Number of controls = 112 Number of selected controls = 102 Wald chi2(3) = 53.15 Prob > chi2 = 0.0000 (Std. err. adjusted for 80 clusters in country_id) ------------------------------------------------------------------------------ | Robust ln_co2 | Coefficient std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- ln_gdp | -7.433319 1.628321 -4.57 0.000 -10.62477 -4.241868 ln_gdp_sq | .8401567 .1713522 4.90 0.000 .5043126 1.176001 ln_gdp_cb | -.0310764 .005952 -5.22 0.000 -.0427421 -.0194107 ------------------------------------------------------------------------------ Note: Chi-squared test is a Wald test of the coefficients of the variables of interest jointly equal to zero. Lassos select controls for model estimation. Type lassoinfo to see number of selected variables in each lasso. Note: Lassos are performed accounting for clusters in country_id. . . display _newline "DSL completed: $S_DATE $S_TIME" DSL completed: 31 Mar 2026 15:39:11 . estimates store dsl_plugin . . *---------------------------------------------* . * 5a. Extract DSL coefficients * . *---------------------------------------------* . matrix dsl_coefs = e(b) . . local b1_dsl = dsl_coefs[1,1] . local b2_dsl = dsl_coefs[1,2] . local b3_dsl = dsl_coefs[1,3] . . display _newline "=== DSL Coefficient Signs ===" === DSL Coefficient Signs === . display " b1 (GDP): " %9.4f `b1_dsl' _col(35) cond(`b1_dsl' < 0, "NEGATIVE > ", "positive") b1 (GDP): -7.4333 NEGATIVE . display " b2 (GDP^2): " %9.4f `b2_dsl' _col(35) cond(`b2_dsl' > 0, "POSITIVE > ", "negative") b2 (GDP^2): 0.8402 POSITIVE . display " b3 (GDP^3): " %9.4f `b3_dsl' _col(35) cond(`b3_dsl' < 0, "NEGATIVE > ", "positive") b3 (GDP^3): -0.0311 NEGATIVE . . local se1_dsl = sqrt(e(V)[1,1]) . local se2_dsl = sqrt(e(V)[2,2]) . local se3_dsl = sqrt(e(V)[3,3]) . . *---------------------------------------------* . * 5b. DSL turning points * . *---------------------------------------------* . local disc_dsl = (`b2_dsl')^2 - 3 * `b1_dsl' * `b3_dsl' . . display _newline "=== DSL Turning Points ===" === DSL Turning Points === . display " Discriminant: " %8.4f `disc_dsl' /// > cond(`disc_dsl' > 0, " (turning points exist)", " (no turning points)") Discriminant: 0.0129 (turning points exist) . . local min_dsl = exp((-`b2_dsl' + sqrt(`disc_dsl')) / (3 * `b3_dsl')) . local max_dsl = exp((-`b2_dsl' - sqrt(`disc_dsl')) / (3 * `b3_dsl')) . . display " Minimum: $" %8.0fc `min_dsl' Minimum: $ 2,429 . display " Maximum: $" %8.0fc `max_dsl' Maximum: $ 27,672 . . *---------------------------------------------* . * 5c. LASSO selection summary * . *---------------------------------------------* . display _newline "=== LASSO Selection ===" === LASSO Selection === . lassoinfo Estimate: active Command: dsregress ------------------------------------------------------ | No. of | Selection selected Variable | Model method lambda variables ------------+----------------------------------------- ln_co2 | linear plugin .3818852 102 ln_gdp | linear plugin .3818852 100 ln_gdp_sq | linear plugin .3818852 100 ln_gdp_cb | linear plugin .3818852 100 ------------------------------------------------------ . . *---------------------------------------------* . * 5d. Pooled DSL (without fixed effects) * . *---------------------------------------------* . * Run DSL without country/year FE to show how LASSO behaves . * when it has only 12 candidate controls (no FE dummies). . * This demonstrates LASSO's selection power but also the cost . * of omitting fixed effects (severe omitted variable bias). . . display _newline "=== DSL WITHOUT FIXED EFFECTS (POOLED) ===" === DSL WITHOUT FIXED EFFECTS (POOLED) === . dsregress $outcome $gdp_vars, /// > controls($controls) /// > vce(cluster country_id) Estimating lasso for ln_co2 using plugin Estimating lasso for ln_gdp using plugin Estimating lasso for ln_gdp_sq using plugin Estimating lasso for ln_gdp_cb using plugin Double-selection linear model Number of obs = 1,600 Number of controls = 12 Number of selected controls = 7 Wald chi2(3) = 25.05 Prob > chi2 = 0.0000 (Std. err. adjusted for 80 clusters in country_id) ------------------------------------------------------------------------------ | Robust ln_co2 | Coefficient std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- ln_gdp | -22.03297 5.277295 -4.18 0.000 -32.37628 -11.68966 ln_gdp_sq | 2.366878 .5652276 4.19 0.000 1.259052 3.474703 ln_gdp_cb | -.084224 .0199055 -4.23 0.000 -.1232381 -.04521 ------------------------------------------------------------------------------ Note: Chi-squared test is a Wald test of the coefficients of the variables of interest jointly equal to zero. Lassos select controls for model estimation. Type lassoinfo to see number of selected variables in each lasso. Note: Lassos are performed accounting for clusters in country_id. . estimates store dsl_pooled . . matrix dsl_pooled_coefs = e(b) . local b1_pooled = dsl_pooled_coefs[1,1] . local b2_pooled = dsl_pooled_coefs[1,2] . local b3_pooled = dsl_pooled_coefs[1,3] . . local se1_pooled = sqrt(e(V)[1,1]) . local se2_pooled = sqrt(e(V)[2,2]) . local se3_pooled = sqrt(e(V)[3,3]) . . display _newline "Pooled DSL coefficients:" Pooled DSL coefficients: . display " b1 = " %9.4f `b1_pooled' b1 = -22.0330 . display " b2 = " %9.4f `b2_pooled' b2 = 2.3669 . display " b3 = " %9.4f `b3_pooled' b3 = -0.0842 . . * Turning points . local disc_pooled = (`b2_pooled')^2 - 3 * `b1_pooled' * `b3_pooled' . if `disc_pooled' > 0 { . local min_pooled = exp((-`b2_pooled' + sqrt(`disc_pooled')) / (3 * `b3_po > oled')) . local max_pooled = exp((-`b2_pooled' - sqrt(`disc_pooled')) / (3 * `b3_po > oled')) . display " Minimum: $" %8.0fc `min_pooled' Minimum: $ 5,581 . display " Maximum: $" %8.0fc `max_pooled' Maximum: $ 24,532 . } . . lassoinfo Estimate: active Command: dsregress ------------------------------------------------------ | No. of | Selection selected Variable | Model method lambda variables ------------+----------------------------------------- ln_co2 | linear plugin .3818852 5 ln_gdp | linear plugin .3818852 7 ln_gdp_sq | linear plugin .3818852 7 ln_gdp_cb | linear plugin .3818852 7 ------------------------------------------------------ . . . *============================================================================ > =* . * SECTION 6: COMPARISON . *============================================================================ > =* . . display _newline "=============================================" ============================================= . display " COMPARISON: ALL METHODS" COMPARISON: ALL METHODS . display "=============================================" ============================================= . display _newline " Sparse FE Kitchen FE BMA (UIP) DSL" Sparse FE Kitchen FE BMA (UIP) DSL . display " b1 (GDP): " %9.4f `b1_sparse' " " %9.4f `b1_kitchen' " " %9.4f ` > b1_bma' " " %9.4f `b1_dsl' b1 (GDP): -7.4980 -7.1307 -7.1390 -7.4333 . display " b2 (GDP^2): " %9.4f `b2_sparse' " " %9.4f `b2_kitchen' " " %9.4f ` > b2_bma' " " %9.4f `b2_dsl' b2 (GDP^2): 0.8490 0.8060 0.8078 0.8402 . display " b3 (GDP^3): " %9.4f `b3_sparse' " " %9.4f `b3_kitchen' " " %9.4f ` > b3_bma' " " %9.4f `b3_dsl' b3 (GDP^3): -0.0315 -0.0298 -0.0299 -0.0311 . . display _newline " TRUE DGP: b1 = -7.1000, b2 = 0.8100, b3 = -0.0300" TRUE DGP: b1 = -7.1000, b2 = 0.8100, b3 = -0.0300 . . *---------------------------------------------* . * 6a. Predicted EKC curves -- Figure 6 * . *---------------------------------------------* . * Normalize both curves so they are comparable on the same y-scale . . quietly summarize ln_gdp . local xmin = r(min) . local xmax = r(max) . local xmean = r(mean) . . preserve . clear . set obs 500 Number of observations (_N) was 0, now 500. . gen lngdp = `xmin' + (_n - 1) * (`xmax' - `xmin') / 499 . . * Compute cubic component for each method . gen fit_bma = `b1_bma' * lngdp + `b2_bma' * lngdp^2 + `b3_bma' * lngdp^3 . gen fit_dsl = `b1_dsl' * lngdp + `b2_dsl' * lngdp^2 + `b3_dsl' * lngdp^3 . . * Normalize: subtract value at the sample mean GDP . local norm_bma = `b1_bma' * `xmean' + `b2_bma' * `xmean'^2 + `b3_bma' * `xmea > n'^3 . local norm_dsl = `b1_dsl' * `xmean' + `b2_dsl' * `xmean'^2 + `b3_dsl' * `xmea > n'^3 . . replace fit_bma = fit_bma - `norm_bma' (500 real changes made) . replace fit_dsl = fit_dsl - `norm_dsl' (500 real changes made) . . * Turning point vertical lines (in log scale) . local lnmin_bma = ln(`min_bma') . local lnmax_bma = ln(`max_bma') . local lnmin_dsl = ln(`min_dsl') . local lnmax_dsl = ln(`max_dsl') . . twoway /// > (line fit_bma lngdp, lcolor("106 155 204") lwidth(medthick) lpattern(soli > d)) /// > (line fit_dsl lngdp, lcolor("217 119 87") lwidth(medthick) lpattern(dash) > ), /// > xline(`lnmin_bma', lcolor("106 155 204"%50) lpattern(shortdash) lwidth(th > in)) /// > xline(`lnmax_bma', lcolor("106 155 204"%50) lpattern(shortdash) lwidth(th > in)) /// > xline(`lnmin_dsl', lcolor("217 119 87"%50) lpattern(shortdash) lwidth(thi > n)) /// > xline(`lnmax_dsl', lcolor("217 119 87"%50) lpattern(shortdash) lwidth(thi > n)) /// > yline(0, lcolor(gs12) lpattern(dot)) /// > ytitle("Predicted log CO2 (normalized at mean GDP)") /// > xtitle("Log GDP per capita") /// > title("Predicted EKC Shape: BMA vs. DSL", size(medium)) /// > subtitle("Both methods trace an inverted-N curve", size(small)) /// > legend(order(1 "BMA" 2 "DSL") rows(1) position(6) size(small)) /// > note("Curves normalized to zero at sample-mean GDP." /// > "Vertical dashed lines mark turning points (blue = BMA, orange = DSL > ).") /// > scheme(s2color) /// > name(fig5_ekc, replace) . . graph export "stata_bma_dsl_fig5_ekc_curves.png", replace width(2400) file stata_bma_dsl_fig5_ekc_curves.png written in PNG format . display _newline "Saved: stata_bma_dsl_fig5_ekc_curves.png" Saved: stata_bma_dsl_fig5_ekc_curves.png . restore . . *---------------------------------------------* . * 6b. Answer-key evaluation -- Figure 7 * . *---------------------------------------------* . * Compare BMA PIPs against ground truth . . estimates restore bma_uip (results bma_uip are active now) . matrix pip_mat = e(pip) . local nvars_all = colsof(pip_mat) . . preserve . clear . set obs `nvars_all' Number of observations (_N) was 0, now 116. . . gen varname = "" (116 missing values generated) . gen pip = . (116 missing values generated) . . local varnames : colnames pip_mat . forvalues i = 1/`nvars_all' { 2. local vname : word `i' of `varnames' 3. replace varname = "`vname'" in `i' 4. replace pip = pip_mat[1,`i'] in `i' 5. } variable varname was str1 now str6 (1 real change made) (1 real change made) variable varname was str6 now str9 (1 real change made) (1 real change made) (1 real change made) (1 real change made) variable varname was str9 now str11 (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) variable varname was str11 now str13 (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) (1 real change made) . . * Keep ONLY the 15 candidate variables . keep if inlist(varname, "ln_gdp", "ln_gdp_sq", "ln_gdp_cb", /// > "fossil_fuel", "renewable", "urban", "globalization") | /// > inlist(varname, "pop_density", "democracy", "corruption", /// > "industry", "services", "trade", "fdi", "credit") (101 observations deleted) . local nvars = _N . . * Readable labels . gen label = "" (15 missing values generated) . replace label = "GDP per capita (log)" if varname == "ln_gdp" variable label was str1 now str20 (1 real change made) . replace label = "GDP squared (log)" if varname == "ln_gdp_sq" (1 real change made) . replace label = "GDP cubed (log)" if varname == "ln_gdp_cb" (1 real change made) . replace label = "Fossil fuel share (%)" if varname == "fossil_fuel" variable label was str20 now str21 (1 real change made) . replace label = "Renewable energy (%)" if varname == "renewable" (1 real change made) . replace label = "Urban population (%)" if varname == "urban" (1 real change made) . replace label = "Globalization index" if varname == "globalization" (1 real change made) . replace label = "Population density" if varname == "pop_density" (1 real change made) . replace label = "Democracy score" if varname == "democracy" (1 real change made) . replace label = "Corruption index" if varname == "corruption" (1 real change made) . replace label = "Industry VA (% GDP)" if varname == "industry" (1 real change made) . replace label = "Services VA (% GDP)" if varname == "services" (1 real change made) . replace label = "Trade openness (% GDP)" if varname == "trade" variable label was str21 now str22 (1 real change made) . replace label = "FDI inflows (% GDP)" if varname == "fdi" (1 real change made) . replace label = "Domestic credit (% GDP)" if varname == "credit" variable label was str22 now str23 (1 real change made) . . * Ground truth . gen is_true = inlist(varname, "fossil_fuel", "renewable", "urban", /// > "democracy", "industry", "ln_gdp", "ln_gdp_sq", "ln_gdp_cb") . . * Sort: true predictors first (by PIP), then noise (by PIP) . gsort -is_true -pip . gen order = _n . labmask order, values(label) . . * Dot plot: BMA PIP by variable, colored by ground truth . graph twoway /// > (scatter order pip if is_true == 1, /// > mcolor("106 155 204") msymbol(circle) msize(large)) /// > (scatter order pip if is_true == 0, /// > mcolor(gs9) msymbol(diamond) msize(large)), /// > xline(0.8, lcolor("217 119 87") lpattern(dash) lwidth(medium)) /// > ylabel(1(1)`nvars', valuelabel angle(0) labsize(small) nogrid) /// > xlabel(0(0.2)1, format(%3.1f)) /// > ytitle("") /// > xtitle("BMA Posterior Inclusion Probability") /// > title("Answer Key: Do BMA and DSL Recover the Truth?", size(medsmall)) // > / > subtitle("True predictors should have PIP > 0.8; noise should have PIP < > 0.8", size(small)) /// > legend(order(1 "True predictor" 2 "Noise variable") /// > rows(1) position(6) size(vsmall)) /// > note("Dashed line = 0.8 threshold. Circle = true predictor, diamond = noi > se." /// > "Only candidate variables shown (country and year FE excluded).", si > ze(vsmall)) /// > scheme(s2color) ysize(7) xsize(9) /// > name(fig6_answer, replace) . . graph export "stata_bma_dsl_fig6_answer_key.png", replace width(2400) file stata_bma_dsl_fig6_answer_key.png written in PNG format . display _newline "Saved: stata_bma_dsl_fig6_answer_key.png" Saved: stata_bma_dsl_fig6_answer_key.png . restore . . *---------------------------------------------* . * 6c. Export comparison CSV * . *---------------------------------------------* . preserve . clear . set obs 4 Number of observations (_N) was 0, now 4. . . gen method = "" (4 missing values generated) . replace method = "Sparse FE" in 1 variable method was str1 now str9 (1 real change made) . replace method = "Kitchen-Sink FE" in 2 variable method was str9 now str15 (1 real change made) . replace method = "BMA (UIP)" in 3 (1 real change made) . replace method = "DSL (Plugin)" in 4 (1 real change made) . . gen b_gdp = . (4 missing values generated) . replace b_gdp = `b1_sparse' in 1 (1 real change made) . replace b_gdp = `b1_kitchen' in 2 (1 real change made) . replace b_gdp = `b1_bma' in 3 (1 real change made) . replace b_gdp = `b1_dsl' in 4 (1 real change made) . . gen b_gdp_sq = . (4 missing values generated) . replace b_gdp_sq = `b2_sparse' in 1 (1 real change made) . replace b_gdp_sq = `b2_kitchen' in 2 (1 real change made) . replace b_gdp_sq = `b2_bma' in 3 (1 real change made) . replace b_gdp_sq = `b2_dsl' in 4 (1 real change made) . . gen b_gdp_cb = . (4 missing values generated) . replace b_gdp_cb = `b3_sparse' in 1 (1 real change made) . replace b_gdp_cb = `b3_kitchen' in 2 (1 real change made) . replace b_gdp_cb = `b3_bma' in 3 (1 real change made) . replace b_gdp_cb = `b3_dsl' in 4 (1 real change made) . . gen min_tp = . (4 missing values generated) . replace min_tp = `min_sparse' in 1 (1 real change made) . replace min_tp = `min_kitchen' in 2 (1 real change made) . replace min_tp = `min_bma' in 3 (1 real change made) . replace min_tp = `min_dsl' in 4 (1 real change made) . . gen max_tp = . (4 missing values generated) . replace max_tp = `max_sparse' in 1 (1 real change made) . replace max_tp = `max_kitchen' in 2 (1 real change made) . replace max_tp = `max_bma' in 3 (1 real change made) . replace max_tp = `max_dsl' in 4 (1 real change made) . . export delimited "stata_bma_dsl_comparison.csv", replace file stata_bma_dsl_comparison.csv saved . display _newline "Saved: stata_bma_dsl_comparison.csv" Saved: stata_bma_dsl_comparison.csv . restore . . . *============================================================================ > =* . * APPENDIX A: FIRST-DIFFERENCES ANALYSIS . *============================================================================ > =* . * Create a cross-sectional dataset by taking (2014 value) - (1995 value) . * for each country. This removes time-invariant FE and produces a setting . * where BMA and DSL operate on pure cross-sectional data (N=80). . . display _newline "=============================================" ============================================= . display " APPENDIX A: FIRST DIFFERENCES" APPENDIX A: FIRST DIFFERENCES . display "=============================================" ============================================= . . preserve . . * Keep only first and last years . keep if year == 1995 | year == 2014 (1,440 observations deleted) . . * Reshape to wide . reshape wide $outcome $gdp_vars $controls, i(country_id) j(year) (j = 1995 2014) Data Long -> Wide ----------------------------------------------------------------------------- Number of observations 160 -> 80 Number of variables 23 -> 38 j variable (2 values) year -> (dropped) xij variables: ln_co2 -> ln_co21995 ln_co22014 ln_gdp -> ln_gdp1995 ln_gdp2014 ln_gdp_sq -> ln_gdp_sq1995 ln_gdp_sq2014 ln_gdp_cb -> ln_gdp_cb1995 ln_gdp_cb2014 fossil_fuel -> fossil_fuel1995 fossil_fuel2014 renewable -> renewable1995 renewable2014 urban -> urban1995 urban2014 globalization -> globalization1995 globalization2 > 014 pop_density -> pop_density1995 pop_density2014 democracy -> democracy1995 democracy2014 corruption -> corruption1995 corruption2014 industry -> industry1995 industry2014 services -> services1995 services2014 trade -> trade1995 trade2014 fdi -> fdi1995 fdi2014 credit -> credit1995 credit2014 ----------------------------------------------------------------------------- . . * Compute first differences: delta_var = var(2014) - var(1995) . foreach v in $outcome $gdp_vars $controls { 2. gen d_`v' = `v'2014 - `v'1995 3. } . . summarize d_* Variable | Obs Mean Std. dev. Min Max -------------+--------------------------------------------------------- d_ln_co2 | 80 -.2472333 .2530391 -.9905968 .2878819 d_ln_gdp | 80 .3805752 .1730117 .0108681 .7278786 d_ln_gdp_sq | 80 7.280438 3.475221 .1978302 15.87447 d_ln_gdp_cb | 80 106.4234 57.24153 2.700623 259.7539 d_fossil_f~l | 80 -4.725521 4.096132 -16.34528 4.306183 -------------+--------------------------------------------------------- d_renewable | 80 7.752427 2.781472 .6046658 14.84042 d_urban | 80 5.410538 2.27026 -.3689499 10.54026 d_globaliz~n | 80 4.101481 4.485683 -9.696026 14.74178 d_pop_dens~y | 80 22.64725 38.11427 -8.772725 269.2358 d_democracy | 80 -.075163 .8260154 -2.001941 1.812005 -------------+--------------------------------------------------------- d_corruption | 80 -.6094901 6.946694 -12.77691 20.62741 d_industry | 80 -1.976088 2.89653 -10.21347 4.43928 d_services | 80 3.76396 2.84419 -2.819748 11.15166 d_trade | 80 .169289 7.233232 -19.90265 21.00221 d_fdi | 80 .682368 3.022733 -5.047574 9.474737 -------------+--------------------------------------------------------- d_credit | 80 5.748046 7.335664 -9.391449 26.14455 . . *---------------------------------------------* . * A1. FD: Sparse OLS * . *---------------------------------------------* . display _newline "=== FD: Sparse OLS ===" === FD: Sparse OLS === . regress d_ln_co2 d_ln_gdp d_ln_gdp_sq d_ln_gdp_cb, robust Linear regression Number of obs = 80 F(3, 76) = 6.13 Prob > F = 0.0009 R-squared = 0.1433 Root MSE = .23879 ------------------------------------------------------------------------------ | Robust d_ln_co2 | Coefficient std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- d_ln_gdp | -10.36189 4.092422 -2.53 0.013 -18.51265 -2.211121 d_ln_gdp_sq | 1.155962 .4223643 2.74 0.008 .3147506 1.997173 d_ln_gdp_cb | -.0414947 .0143721 -2.89 0.005 -.0701192 -.0128702 _cons | -.3036562 .0724366 -4.19 0.000 -.4479262 -.1593861 ------------------------------------------------------------------------------ . . *---------------------------------------------* . * A2. FD: Kitchen-sink OLS * . *---------------------------------------------* . display _newline "=== FD: Kitchen-sink OLS ===" === FD: Kitchen-sink OLS === . regress d_ln_co2 d_ln_gdp d_ln_gdp_sq d_ln_gdp_cb /// > d_fossil_fuel d_renewable d_urban d_industry d_democracy /// > d_services d_trade d_fdi d_credit d_pop_density /// > d_corruption d_globalization, robust Linear regression Number of obs = 80 F(15, 64) = 2.71 Prob > F = 0.0029 R-squared = 0.3707 Root MSE = .22301 ------------------------------------------------------------------------------ | Robust d_ln_co2 | Coefficient std. err. t P>|t| [95% conf. interval] -------------+---------------------------------------------------------------- d_ln_gdp | -8.109709 5.031758 -1.61 0.112 -18.1618 1.942382 d_ln_gdp_sq | .9238864 .5213262 1.77 0.081 -.1175823 1.965355 d_ln_gdp_cb | -.0336221 .0179583 -1.87 0.066 -.0694979 .0022536 d_fossil_f~l | .0147108 .0067313 2.19 0.033 .0012635 .0281582 d_renewable | -.0237808 .0110384 -2.15 0.035 -.0458327 -.001729 d_urban | .0002501 .014913 0.02 0.987 -.0295421 .0300424 d_industry | .0309085 .0105974 2.92 0.005 .0097377 .0520793 d_democracy | .019337 .0290345 0.67 0.508 -.038666 .07734 d_services | -.0047239 .0098816 -0.48 0.634 -.0244647 .0150169 d_trade | .006726 .0044062 1.53 0.132 -.0020764 .0155284 d_fdi | .0000124 .0091898 0.00 0.999 -.0183463 .0183712 d_credit | .0028644 .0043456 0.66 0.512 -.0058169 .0115457 d_pop_dens~y | .0006396 .0004991 1.28 0.205 -.0003575 .0016366 d_corruption | -.0036115 .0033497 -1.08 0.285 -.0103033 .0030803 d_globaliz~n | -.0004567 .0082494 -0.06 0.956 -.0169368 .0160235 _cons | -.0085823 .1746184 -0.05 0.961 -.3574226 .340258 ------------------------------------------------------------------------------ . . *---------------------------------------------* . * A3. FD: BMA * . *---------------------------------------------* . display _newline "=== FD: BMA ===" === FD: BMA === . bmaregress d_ln_co2 d_ln_gdp d_ln_gdp_sq d_ln_gdp_cb /// > d_fossil_fuel d_renewable d_urban d_industry d_democracy /// > d_services d_trade d_fdi d_credit d_pop_density /// > d_corruption d_globalization, /// > mprior(uniform) gprior(uip) /// > mcmcsize(50000) rseed(9988) pipcutoff(0.5) burnin(5000) Burn-in ... Simulation ... Computing model probabilities ... Bayesian model averaging No. of obs = 80 Linear regression No. of predictors = 15 MC3 sampling Groups = 15 Always = 0 No. of models = 2,317 For CPMP >= .9 = 581 Priors: Mean model size = 3.304 Models: Uniform Burn-in = 5,000 Cons.: Noninformative MCMC sample size = 50,000 Coef.: Zellner's g Acceptance rate = 0.3080 g: Unit-information, g = 80 Shrinkage, g/(1+g) = 0.9877 sigma2: Noninformative Mean sigma2 = 0.051 Sampling correlation = 0.9958 ------------------------------------------------------------------------------ d_ln_co2 | Mean Std. dev. Group PIP -------------+---------------------------------------------------------------- d_industry | .0364834 .0090778 7 .99823 -------------+---------------------------------------------------------------- Always | _cons | -.1551343 .099996 0 1 ------------------------------------------------------------------------------ Note: Coefficient posterior means and std. dev. estimated from 2,317 models. Note: 14 predictors with PIP less than .5 not shown. . . * List all PIPs . matrix pip_fd = e(pip) . local varnames_fd : colnames pip_fd . local ncols_fd = colsof(pip_fd) . display _newline "FD BMA PIPs:" FD BMA PIPs: . forvalues i = 1/`ncols_fd' { 2. local vname : word `i' of `varnames_fd' 3. display " `vname': " %8.6f pip_fd[1,`i'] 4. } d_ln_gdp: 0.298283 d_ln_gdp_sq: 0.267234 d_ln_gdp_cb: 0.271081 d_fossil_fuel: 0.183373 d_renewable: 0.349742 d_urban: 0.096227 d_industry: 0.998225 d_democracy: 0.094103 d_services: 0.095686 d_trade: 0.179387 d_fdi: 0.089272 d_credit: 0.092881 d_pop_density: 0.100171 d_corruption: 0.099431 d_globalization: 0.089335 _cons: 1.000000 . . *---------------------------------------------* . * A4. FD: DSL * . *---------------------------------------------* . display _newline "=== FD: DSL ===" === FD: DSL === . dsregress d_ln_co2 d_ln_gdp d_ln_gdp_sq d_ln_gdp_cb, /// > controls(d_fossil_fuel d_renewable d_urban d_industry d_democracy /// > d_services d_trade d_fdi d_credit d_pop_density /// > d_corruption d_globalization) /// > rseed(9988) Estimating lasso for d_ln_co2 using plugin Estimating lasso for d_ln_gdp using plugin Estimating lasso for d_ln_gdp_sq using plugin Estimating lasso for d_ln_gdp_cb using plugin Double-selection linear model Number of obs = 80 Number of controls = 12 Number of selected controls = 1 Wald chi2(3) = 10.65 Prob > chi2 = 0.0138 ------------------------------------------------------------------------------ | Robust d_ln_co2 | Coefficient std. err. z P>|z| [95% conf. interval] -------------+---------------------------------------------------------------- d_ln_gdp | -5.047196 4.558593 -1.11 0.268 -13.98187 3.887483 d_ln_gdp_sq | .5943786 .4700569 1.26 0.206 -.326916 1.515673 d_ln_gdp_cb | -.0220809 .0160386 -1.38 0.169 -.0535159 .0093541 ------------------------------------------------------------------------------ Note: Chi-squared test is a Wald test of the coefficients of the variables of interest jointly equal to zero. Lassos select controls for model estimation. Type lassoinfo to see number of selected variables in each lasso. . . lassoinfo Estimate: active Command: dsregress ------------------------------------------------------ | No. of | Selection selected Variable | Model method lambda variables ------------+----------------------------------------- d_ln_co2 | linear plugin .3818852 1 d_ln_gdp | linear plugin .3818852 0 d_ln_gdp_sq | linear plugin .3818852 0 d_ln_gdp_cb | linear plugin .3818852 0 ------------------------------------------------------ . . restore . . . *============================================================================ > =* . * WRAP-UP . *============================================================================ > =* . . display _newline "=============================================" ============================================= . display " TUTORIAL COMPLETE: $S_DATE $S_TIME" TUTORIAL COMPLETE: 31 Mar 2026 15:39:14 . display "=============================================" ============================================= . display _newline "Output files:" Output files: . display " analysis.log -- this log" analysis.log -- this log . display " stata_bma_dsl_fig1_scatter.png -- scatter plot" stata_bma_dsl_fig1_scatter.png -- scatter plot . display " stata_bma_dsl_fig2_instability.png -- coefficient instability" stata_bma_dsl_fig2_instability.png -- coefficient instability . display " stata_bma_dsl_fig3_pip.png -- BMA PIPs (color-coded)" stata_bma_dsl_fig3_pip.png -- BMA PIPs (color-coded) . display " stata_bma_dsl_fig4_coefdensity.png -- BMA coefficient densities (4 > panels)" stata_bma_dsl_fig4_coefdensity.png -- BMA coefficient densities (4 panels) . display " stata_bma_dsl_fig5_ekc_curves.png -- EKC curves (normalized)" stata_bma_dsl_fig5_ekc_curves.png -- EKC curves (normalized) . display " stata_bma_dsl_fig6_answer_key.png -- answer key evaluation" stata_bma_dsl_fig6_answer_key.png -- answer key evaluation . display " stata_bma_dsl_comparison.csv -- coefficient comparison" stata_bma_dsl_comparison.csv -- coefficient comparison . . capture erase "_bma_temp.dta" . log close name: log: /Users/carlosmendez/Documents/Github/starter-academic-v501/content > /post/stata_bma_dsl/analysis.log log type: text closed on: 31 Mar 2026, 15:39:14 -------------------------------------------------------------------------------