program define iepowerdemo, rclass
*!	version 0.5 4Feb2026 by Mead Over for UPSE DE 292
    version 14
    syntax , Reps(integer) n1(integer) n2(integer)   ///  Required options
		[mu1(real 0) sigma1(real 1)   ///
		 mu2(real 1) sigma2(real 1)   ///   Optional options
		title(passthru) xtitle(passthru)  ///
		level(integer 95) Filestub(string) Verbose SHOWrep(integer 0)]

	cap which require
	if _rc==0 {
		require hist_overlay>=1.41, install from(http://digital.cgdev.org/doc/stata/MO/Misc)
	}

	if "`verbose'"~=""  {
		local shhh noisily
	}
	else {
		local shhh quietly
	}

	local Ntotal = `n1' + `n2'

	if "`filestub'"=="" {
		local filestub obs_`Ntotal'
	}

	cap assert `showrep' >= 0 & `showrep' <= `reps'
	if _rc>0 {
		di as err "Option -showrep(integer)- must be between 0 and `reps'"
		exit 198
	}

    // Create a holding place named -tempspace- to store the means,
	// their difference and other results from each replication.
	// Specify the DTA file to be used to store the -tempspace- data
	// to the hard disk.  Give the DTA file a temporary filename
	// to assure that it disappears when this ADO file terminates.
    tempname tempspace
	tempfile tmpfilename

	// Open a holding place in which each observation contains the results of one replication.
    `shhh' postfile `tempspace' long(rep Ntotal)  ///
	    long(N_Trt)    float(Mean_Trt SD_Trt)     ///  <-- Treatment Group
	    long(N_Ctrl)   float(Mean_Ctrl SD_Ctrl)   ///  <-- Control Group
	    float(Diff SE_Diff) long(df_t)            ///  Diff in Means, SE of Dif, DF
	    float(t_stat) float(p_value) byte(reject) ///  test results
		using `tmpfilename', replace

    // Run the simulation
    forvalues i = 1/`reps' {
		if "`verbose'"=="" & `i'==`showrep' {
			local shhh noisily
		}
		`shhh' di _n as txt "******************* Replication " as res "`i'" as txt " *******************"
        quietly {

            // Generate random sample with mean=Mu1 and SD=sigma1
            drop _all
            set obs `Ntotal'
			gen long index = _n
				lab var index "Sequence number of observation"

			// Assign the first n1 observations to the Control Group
			// and the rest to the treatment group
			gen byte grp = cond(index<=`n1', 0, 1)
				label values grp grp

			gen y = rnormal(`mu1', `sigma1') if grp ==0
				replace y = rnormal(`mu2', `sigma2') if grp ==1

			// Test the hypothesis that the Control & Treatment means are different

			//   Stata's -ttest- command does all the work. It computes the mean
			//   and standard deviation in each group, the difference between the means
			//   in the two groups, the standard error of the difference in means
			//   and the t-statistic of the test that the test of the null hypothesis
			//   that the difference is zero.

			//   After making one small adjustment, all we need to do is to post
			//   the results from -ttest- to the -tempspace- for each replication.

			// 	 By default, Stata's -ttest- command defines the difference in means
			// 	 as the mean of group 1 minus the mean of group 2.  In contrast,
			// 	 the impact evaluation convention is to define group 1 as the control group
			// 	 and group 2 as the treatment group and the impact as the difference
			// 	 calculated as the treatment mean minus the control mean.
			// 	 To convert -ttest-'s results to the IE convention, we use option -reverse-
			// 	 on the -ttest- command and then rename the results in the posted file.

			if "`verbose'"=="verbose" | `i'==`showrep' {
				noi di _n as txt "Simulated data for replication " as res "`i'"  ///
					_n as inp "summarize"
				noi summarize
				noi di _n as inp "ttest y , by(grp) unequal reverse"
			}

			`shhh' {
				ttest y , by(grp) unequal reverse
					local Trt_Minus_Cntrl = `r(mu_1)' -  `r(mu_2)'
				// return list
			}

			local signif = (100 - `level')/100

			local reject = cond(`r(p)'<`signif',1,0)

			// Post the results from the i'th replication
            post `tempspace' (`i') (`Ntotal')   ///
				(`r(N_1)') (`r(mu_1)') (`r(sd_1)')         ///  <-- Treatment Group
				(`r(N_2)') (`r(mu_2)') (`r(sd_2)')         ///  <-- Control Group
				(`Trt_Minus_Cntrl') (`r(se)') (`r(df_t)')  ///  Diff in Means, SE of Diff, df of t
				(`r(t)') (`r(p)') (`reject')               //   t-stat, p-value

        }
		if "`verbose'"=="" & `i'==`showrep' {
			local shhh quietly
		}
    }  //  End of the simulation loop over -reps- replications

	// Close the "postfile"
    `shhh' postclose `tempspace'

    // Annotate the posted data and save it to the file -sample_means_wide_`filestub'.dta-
    use `tmpfilename', clear

	lab var rep       "Replication sequence number"
	lab var Ntotal    "Total sample size"
	lab var N_Trt     "Sample size for treatment group"
	lab var Mean_Trt  "Mean outcome for treatment group"
	lab var SD_Trt    "SD of outcome for treatment group"
	lab var N_Ctrl    "Sample size for control group"
	lab var Mean_Ctrl "Mean outcome for control group"
	lab var SD_Ctrl   "SD of outcome for control group"
	lab var Diff      "Estimated impact = Mean_Trt - Mean_Ctrl"
	lab var SE_Diff   "Standard error of Diff"
	lab var df_t      "Degrees of freedom for test of difference in means"
	lab var t_stat    "t-statistic = Diff/SE_Diff"
	lab var p_value   "p-value for t_stat with DF = df_t"
	lab var reject    "= 1 (Ho of equal means rejected) if p <. `signif'"

	label define reject 0 "Accept Ho (p>=`signif')" 1 "Reject Ho (p<`signif')"
	label values reject reject


	`shhh' sum reject
		local reject_nmbr = `r(sum)'
		local reject_prop = `r(mean)'
		local reject_prcnt = 100*`r(mean)'
		local reject_prop_txt : di %6.3f `reject_prop'
		local reject_prcnt_txt : di %5.1f `reject_prcnt'

	label data "Test of the difference between means over `reps' reps. one obs. per rep."

	// Average estimated impact across the replications
	`shhh' summarize Diff
		local Diff_avg = `r(mean)'
		local Diff_avg_txt : di %6.2f `Diff_avg'


	local TrueDiff = `mu2' - `mu1'
	notes _dta: The true impact is defined as the difference between the two population means: `mu2' - `mu1' = `TrueDiff'
	notes _dta: Sample sizes: Total `Ntotal'; Control group: `n1';  Treatment group: `n2'
	notes _dta: Population Means: Control group: `mu1'.;  Treatment group: `mu2';
	notes _dta: Population SDs: Control group: `sigma1';  Treatment group: `sigma2'
	notes _dta: Average difference between the means: `Diff_avg'
	notes _dta: Out of `reps' replications, at a dignificance level of alpha = `signif', Ho is rejected `reject_nmbr' times, for a rejection rate of `reject_prcnt_txt'%

	qui compress
    `shhh' save sample_means_wide_`filestub', replace

	// To visualize the sampling distribution of the difference between the means
	// of the two groups, we use Stata's -histogram-
	`shhh' histogram Diff, frequency  ///
        subtitle(`"At a sample size of `Ntotal', Ho is rejected `reject_nmbr' times out of `reps' reps."'  ///
			`"for a rejection rate of `reject_prcnt_txt'%"')  ///
		ylabel(,angle(horizontal))  ///
		xline(`TrueDiff')  ///
        xlabel(,grid) ylabel(,grid) xlabel(#7 `TrueDiff')  ///
		name(diffmeans_`filestub', replace)  saving(diffmeans_`filestub'.gph, replace)

	// To distinguish observations on which the null hypothesis is rejected,
	// we use community-contributed -hist_overlay-.
	if `reject_prcnt_txt' > 0 & `reject_prcnt_txt' < 100 {
		`shhh' hist_overlay Diff, frequency over(reject) ///
			subtitle(`"At a sample size of `Ntotal', Ho is rejected `reject_nmbr' times out of `reps' reps."'  ///
				`"for a rejection rate of `reject_prcnt_txt'%"')  ///
			ylabel(,angle(horizontal))  ///
			color1(red%30) color2(green%30)  ///
			xline(`TrueDiff')  ///
			xlabel(,grid) ylabel(,grid) xlabel(#7 `TrueDiff') ///
			name(reject_ho_`filestub', replace) saving(reject_ho_`filestub'.gph, replace)
	}
	else {
		di as txt "The -hist_overlay- visualization of the Ho rejection rate"  ///
			_n as txt "requires a rate greater than 0% and less than 100%."
		di  as txt "At a sample size of " as res "`Ntotal'"  ///
			as txt ", Ho is rejected " as res "`reject_nmbr'"   ///
			as txt " times out of " as res "`reps'" as txt " replications."  ///
			_n as txt `"for a rejection rate of "' as res "`reject_prcnt_txt'%"  ///
			_n as txt "Therefore, the graph file " as res "`reject_ho_`filestub'.gph'" as txt " is not produced."
	}

	// To visualize the overlapping sampling distributions of the means
	// for the two groups, we use community-contributed -hist_overlay-.
	// This program requires a distinct observation for each distribution.
	// So we reshape the data such that there are two observations
	// for each replication, one for each group mean.

	//  Set up for reshaping the data, so that -hist_overlay- can display
	//  the sampling distributions of the two means.
	if "`verbose'"=="verbose" {
		noi di _n as inp "reshape long N_ Mean_ SD_, i(N_* Mean_* SD_*) j(Group) string"
	}
	reshape clear
	reshape xij N_ SD_ Mean_
    reshape j   Group Ctrl Trt, string
	reshape i   N_* Mean_* SD_*
	reshape xi rep Ntotal Diff SE_Diff t_stat df_t p_value reject

	//  Reshape to two observations for each replication, one for each group
	`shhh' reshape long

	encode Group, gen(grp)
		`shhh' replace grp = grp-1
		recast byte grp, force
		cap label drop grp
		label define grp 0 "Control" 1 "Treatment"
		label values grp grp

	keep rep grp Group Ntotal N_ SD_ Mean_ Diff SE_Diff t_stat df_t p_value reject
	sort rep grp
	order rep grp Group N* Mean* SD*

	lab var Group  "Control or Treatment Group"
	lab var grp    "0 = Control; 1 = Treatment"
	lab var N_     "Sample size on this replication by group"
	lab var Mean_  "Mean outcome on this replication by group"
	lab var SD_    "SD of outcome on this replication by group"

	`shhh' compress
	label data "Test of the difference between means over `reps' reps., two obs. per rep."
	`shhh' save sample_means_long_`filestub', replace

	//  Temporarily redefine the variable name for Mean_
	//  to serve as a label for the x-axis of the histogram.
	lab var Mean_ "Sample means by group in `reps' replications"
	`shhh' hist_overlay Mean_, over(grp) ///
		`title' `xtitle' ///
        subtitle(`"At a sample size of `Ntotal', Ho is rejected `reject_nmbr' times out of `reps' reps."'  ///
			`"for a rejection rate of `reject_prcnt_txt'%"')  ///
		ylabel(,angle(horizontal))  ///
		xline1(`mu1') xline2(`mu2')  ///
        xlabel(,grid) ylabel(,grid) xlabel(#7 `mu1' `mu2')  ///
		name(twomeans_`filestub', replace) saving(twomeans_`filestub'.gph, replace)

	di _n as txt "iepowerdemo has saved the following files to disk:"
	dir *_`filestub'.gph
	dir *_`filestub'.dta
	local gphfiles  : dir . files  "*_`filestub'.gph"
	local datafiles : dir . files  "*_`filestub'.dta"

	return local datafiles `datafiles'
	return local gphfiles  `gphfiles'

end
exit
*	ver 0.1  Make the graph title optional
*	ver 0.2  Add the option -xtitle-
*	ver 0.3  Compute the means of both distributions within the sample from a single replication
*	ver 0.5  4Feb2026 rename program to -iepowerdemo.ado-.  Add option -showrep()-
*				Deprecate the option -xlabel-, in favor of graph specific -xlabel- specifications.