*********************************************************
clear

*** Set sample size to be drawn, means and standard deviations for the two samples
*** to be sampled, and the number of simulation replications
global SampSize 50
global MeanS1 50
global MeanS2 50
global StdDevS1 1
global StdDevS2 1
local RepNum 100

***  NO NEED TO CHANGE ANYTHING BELOW UNLESS YOU KNWO WHAT YOU ARE DOING  ***
matrix m = ($MeanS1, $MeanS2)
matrix sd = ($StdDevS1,$StdDevS2)

*** Define program simttest
capture program drop simttest
program define simttest, rclass
	quietly {
		drawnorm Sa1 Sb1, n($SampSize) means(m) sds(sd)
			*** Summary statistics 
			summarize Sa1
				return scalar mean1=r(mean)
				return scalar sd1=r(sd)
			summarize Sb1
				return scalar mean2=r(mean)
				return scalar sd2=r(sd)
			*** t-test
			ttest Sa1 == Sb1, unpaired
				return scalar t = r(t)
				return scalar p = r(p)
	}
	drop Sa1 Sb1		
end

*** Start simulation ***
simulate "simttest" mean1=r(mean1) sd1=r(sd1) mean2=r(mean2) sd2=r(sd2) ///
	t=r(t) pvalue=r(p), reps(`RepNum') dots
label variable mean1 "Mean of Sample 1"
label variable mean2 "Mean of Sample 2"
label variable sd1 "Standard Deviation of Sample 1"
label variable sd2 "Standard Deviation of Sample 2"

scatter mean2 mean1 if (pvalue > 0.05) || ///
	scatter mean2 mean1 if (pvalue <= 0.05), mcolor(red) ///
	ysize(4.4) xsize(4) ylabel(,grid) xlabel(,grid) ///
	xline($MeanS1, lstyle(foreground) lwidth(med)) ///
	yline($MeanS2, lstyle(foreground) lwidth(med)) ///
	legend(label(1 "p-value > 0.05") label(2 "p-value <= 0.05")) ///
	title(Comparison of two samples using t-test) ///
	subtitle(p-value: Sample 1 Mean == Sample 2 Mean) name(TwoSamples, replace)
scatter sd2 sd1 if (pvalue > 0.05) || ///
	scatter sd2 sd1 if (pvalue <= 0.05), mcolor(red) ///
	ysize(4.4) xsize(4) ylabel(,grid) xlabel(,grid) ///
	xline($StdDevS1, lstyle(foreground) lwidth(med)) ///
	yline($StdDevS2, lstyle(foreground) lwidth(med)) ///
	legend(label(1 "p-value > 0.05") label(2 "p-value <= 0.05")) ///
	title(Comparison of two samples using t-test) ///
	subtitle(p-value: Sample 1 Mean == Sample 2 Mean) name(TwoSamples2, replace)