* code for displaying table of results largely borrowed from -leanout
* some ideas for significant digits from http://www.stata.com/statalist/archive/2010-01/msg00355.html

* designed only for continuous data - it is user responsibility not to perturb ordinal/categorical/orderd or other types of data
* przetestowac na kazdej komendzie estimation in przykladach z helpa
* oczyscic kod i ladnie poprawic
* obsluge bledow dorobic - jakich bledow?

*! 1.0.0 Michal Brzezinski 11feb2011
program define compsta, rclass
	version 11

	if "`e(cmd)'" == "" {
		error 301
	}
	local command = e(cmdline)
	if "`e(cmdline)'" == "" {
		dis as error "e(cmdline) not found"
		exit
	}
	syntax ,							/// jakie zrobic abbreviations dla opcji??
		PERVar(varlist numeric) 		/// perturbed variables
		RELDig(numlist >=1 integer) 	/// reliable digits
		[Reps(integer 499), 			///
		GRaph							/// ta opcja rysuje wszystkie
		GRSelected(varlist numeric)		///
		NODots							///
		NOHEader						/// for use if only graph is displayed
		NOTABle							///	for use if only graph is displayed
		Level(cilevel)					/// for testing the % inside CI
		*								/// this is for graph_options; * returns `options', see e.g. glcurve
		]
	qui {
		estimates store org
		tempname b V
		mat `b' = e(b)
		mat `V' = e(V)
		local nobs = e(N)
		if "`e(df_r)'" != "" {
			local type t
			local ci_mult = -invttail(e(df_r), 1-(1-`level'/100)/2)
		}
	else {
		local type z
		local ci_mult = invnormal(1-(1-`level'/100)/2)
	}

// perturbing variables
	local varnum: word count `pervar'
	local dignum : word count `reldig'
	if `varnum'!=`dignum' {
			display as error "pervar() and reldig() must have the same number of elements"
			exit 198
	}
// handling other errors - what if pervar and/or reldig empty or
// misnamed or reldig too big, etc.


tempname memhold perturb_b
tempfile master results
save "`master'"
local j = 1
local coefnum=colsof(`b')
local coefnames: colfullnames e(b)
local coefnames ""
forvalues i=1/`coefnum' {
	local coefnames "`coefnames' b`i'"
}

postfile `memhold' double `coefnames' using "`results'", replace


	if "`nodots'" == "" {
		noi dis
		noi _dots 0, title(Perturbation replications) reps(`reps')
	}

	forvalues i = 1/`reps' {
		preserve
		foreach pdig of numlist `reldig' {
			local pervarname: word `j' of `pervar'
			sum `pervarname', meanonly
			local pvmax = r(max)
			*local reldig = `pdig'
			local left = int(log10(abs(`r(max)'))+1)
			local A = 10^(`pdig'-`left')
			replace `pervarname' = (round(`pervarname'*`A')+ ///
			round(-0.5+0.999*runiform(),0.001))/`A'
			local ++j

					}
	*noi list // shows perturbed variables for a replication
	`command'
	restore
	mat `perturb_b' = e(b)
	*noi matlist `perturb_b''
	local posted_coefnames ""
	forvalues k=1/`coefnum' {
		local posted_coefnames "`posted_coefnames' (`perturb_b'[1,`k'])"
	}
	*noi dis `posted_coefnames'
	post `memhold' `posted_coefnames'
	if "`nodots'" == "" {
		nois _dots `i' 0
	}
	local j = 1

	} // end forvalues 1/reps
	postclose `memhold'

use "`results'", clear
estimates restore org
tempname b V
mat `b' = e(b)
mat `V' = e(V)

// collect descriptive. stats
// calculate other statistics
mat perturb_desc = J(`coefnum', 5,.)
local coln: colnames `b'
matrix rownames perturb_desc = `coln'
matrix colnames perturb_desc = mean median sd min max

forvalues i = 1/`coefnum' {

	// descriptives
	sum b`i',d
	scalar b`i'_mean   = r(mean)
	scalar b`i'_diff   = `b'[1,`i'] - b`i'_mean
	scalar b`i'_sd     = r(sd)
	scalar b`i'_min    = r(min)
	scalar b`i'_max    = r(max)
	scalar b`i'_median = r(p50)

	// return matrix with descriptives
	mat perturb_desc[`i', 1] = b`i'_mean
	mat perturb_desc[`i', 2] = b`i'_median
	mat perturb_desc[`i', 3] = b`i'_sd
	mat perturb_desc[`i', 4] = b`i'_min
	mat perturb_desc[`i', 5] = b`i'_max

	// proportion of pert. coefs. within origincal CI
	scalar lb_`i' = `b'[1,`i'] - `ci_mult'*`V'[`i',`i']^0.5
	scalar ub_`i' = `b'[1,`i'] + `ci_mult'*`V'[`i',`i']^0.5
	count if b`i' > lb_`i' & b`i' < ub_`i'
	scalar prop_in_CI_`i' = r(N)/`reps'*100

	// proportion of pert. coefs. that agree to one digit with unperturbed coef.
	local dd = floor(log10(abs(`b'[1,`i'])))
	local ll = `b'[1,`i'] - 0.5*10^`dd'
	local uu = `b'[1,`i'] + 0.5*10^`dd'
	count if b`i' > `ll' & b`i' < `uu'
	scalar prop_agree_`i' = r(N)/`reps'*100
	*noi dis prop_agree_`i'
}
return mat perturb_desc=perturb_desc

// displaying results

if "`nodots'" == "" {
	if mod(`reps',50)!=0 {
			local di2 noi dis
		}
		else {
			local di2 ""
		}
		`di2'
}

if "`noheader'" == "" {
	noi dis _newline as text  "Perturbation test of computational stability" ///
	_col(56) "Number of obs = " as result %7.0g `nobs'
	noi dis as text _col(56) "Replications  = " as result %7.0g `reps'
	noi dis _newline as text _col(7) "command: " as result "`command'"
	noi dis as text _col(2) "perturb vars: " as text "`pervar'"
	noi dis
}

*local eqs: coleq `b'
*dis "`eqs'"
_ms_eq_info, matrix(`b')
local neq = r(k_eq)
forvalues eq = 1/`neq' {
	local k`eq' = r(k`eq')
	local eqname`eq' `"`r(eq`eq')'"'
	dis `k`eq''
}
dis `neq'
local width = 12


if "`notable'" == "" {

	noi display as text "{hline 13}{c TT}{hline 64}"
	*noi display as text _col(14) "{c |}" _col(49)    "Perturbed Coef."
	noi display as text _col(14) "{c |}" _col(18) "Observed" _col(34) "Mean" _col(48) "SD" _col(68) "%" _col(78) "%"
	noi display as text _col(14) "{c |}" _col(21)  "Coef." _col(29) "Perturbed" _col(41) "Perturbed" _col(54) "P-val" _col(64) "in CI" _col(74) "agree"

	local neq = `neq' - 0`e(k_aux)'

	local coln: colnames `b'
	dis "`coln'"

	// Loop over equations
	local vsquish_h `vsquish'
		local i 1
		forvalues eq = 1/`neq' {
			if `"`eqname`eq''"' != `"model"'      &			///
			   `"`eqname`eq''"' != `"_"'          &			///
			   "`e(baselab)'"!="" & "`eqname`eq''"!="`e(baselab)'" & ///
			   `"`eqname`eq''"' != `"`e(depvar)'"' {
				noi display as result abbrev("`eqname`eq''",12) as text "      {c |}           {c |}"
			}
			if "`e(baselab)'"!="" & "`eqname`eq''"=="`e(baselab)'" {
				noi display as text "{hline 13}{c +}{hline 64}"
				noi display as result abbrev("`eqname`eq''",12) _col(14) "{c |}" as text _col(17) "(base outcome)"
				noi display as text "{hline 13}{c +}{hline 64}"
				local i =`k`eq''+1
				local qui2 qui
			}
			else {
			 noi `qui2' display as text "{hline 13}{c +}{hline 64}"

						// Loop over covariates in equations
						// and display
			local vsquish "vsquish"
			forvalues j = 1/`k`eq'' {

				// moze te wszystkie obliczenia poza Display zrobic
				// calculate tstats and p-values for test: uperturb=mean(pert.)

				scalar tstat`i' = (b`i'_mean - `b'[1,`i'])*`reps'^0.5 / b`i'_sd
				scalar pval`i'  = 2*ttail(`reps', abs(tstat`i'))


				local vname: word `i' of `coln'

				*dis "`vname'"
				_ms_parse_parts "`vname'"

				if r(omit) == 1 {
					if "`r(op)'" == "o" {
						noi di as text %12s abbrev("`r(name)'",12) _col(13) "{c |}" as result _col(18) "(omitted)"
						local omit omit
					}
				}

				noi _ms_display, eq(#`eq') el(`j') width(`width') ///
						  nov1bar `first' `vsquish' noomitted

				if r(output) {
					local bi `b'[1,`i']
					noi di as res _col(17) %9.0g `bi'  ///
					_col(29) %9.0g b`i'_mean		   ///
					_col(41) %9.0g b`i'_sd 		   ///
					_col(54) %5.3f pval`i'			   ///
					_col(63) %6.2f prop_in_CI_`i'  	   ///
					_col(73) %6.2f prop_agree_`i'

				local first
				local vsquish `vsquish_h'
				}


				if !r(output) & "`omit'"!="omit" {
					if r(first) {
						local first first
					}
				}
	local ++i
	}

	}
	}

	if 0`e(k_aux)'==0 {
	noi display as text "{hline 13}{c BT}{hline 64}"
	}
	else {
	noi display as text "{hline 13}{c +}{hline 64}"
	}
		// Loop over auxilary parameters
		forvalues aux = 1/0`e(k_aux)' {
			local eq = `neq' + `aux'

			local bi `b'[1,`i']
			local se = sqrt(`V'[`i',`i'])

			noi di as text %11s abbrev("/`eqname`eq''",12) _col(13) "{c |}" _c
			noi di as res _col(17) %9.0g `bi'  ///
			_col(29) %9.0g b`i'_mean		   ///
			_col(41) %9.0g b`i'_sd  		   ///
			_col(58) %9.0g b`i'_min 		   ///
			_col(70) %9.0g b`i'_max

			local ++i
		}

		if 0`e(k_aux)' > 0 {
			noi di as text "{hline 13}{c BT}{hline 64}"

		}
	// footer

	noi dis as text "Note:  P-value for Ho: Diff. = 0, Ha: Diff. != 0"
	noi dis as text `"Note:  "% in CI" refers to CI for observed coef."'
	noi dis as text `"Note:  "% agree" refers to agreement to one significant digit"'
	//definition of agreement to one significant digit.
	//The definition is as follows:
	//Let U be the unperturbed solution and P be a perturbed solution;
	//then, P is said to agree with U if P falls within
	//the interval U +/- five units in the second significant digit of U.
}

if "`graph'" == "graph" {

	tempfile `coefnames'
	noi di
	noi di in text "Wait, drawing graph ..."
	local combine_name ""

	forvalues i=1/`coefnum' {

		local vname: word `i' of `coln'

		local xl = `b'[1,`i']
		twoway kdensity b`i',  xline(`xl', noextend lstyle(refline)) ytitle("") ///
		xtitle("") title(`vname') nodraw saving(`b`i'.gph')
		local combine_name "`combine_name' `"`b`i''"'"
	}

	graph combine `combine_name'  ///
	,title("Monte Carlo distribution of perturbed coefficients") ///
	note("Note: vertical lines represent unperturbed coefficient(s)") `options'

}

if "`grselected'" != "" {

	local varnum_gr: word count `grselected'
	*local varname_gr: word `j' of `grselected'
	local coefnames ""
	forvalues i=1/`varnum_gr' {
		local coefnames "`coefnames' b`i'"
	}

	tempfile `coefnames'
	noi di
	noi di in text "drawing graph ..."
	local combine_name ""

	forvalues i=1/`varnum_gr' {

		local vname: word `i' of `grselected'
		local cfname ""
		local j = 1
		while "`cfname'" != "`vname'" {
			local cfname: word `j' of `coln'
			local ++j
		}
		local --j
		local xl = `b'[1,`j']
		twoway kdensity b`j',  xline(`xl', noextend lstyle(refline)) ytitle("") ///
		xtitle("") title(`vname') nodraw saving(`b`i'.gph')
		*local cn: word `i' of `coefnames'
		local combine_name "`combine_name' `"`b`i''"'"
	}

	graph combine `combine_name'  ///
	,title("Monte Carlo distribution of perturbed coefficients") ///
	note("Note: vertical lines represent unperturbed coefficient(s)") `options'



}

use "`master'", clear
} //end qui
end

exit

---------------------------
Designing table of results for compsta

Perturbation replications (100)
---
....

Perturbation test of computational stability           Number of obs =      74
												       Replications  =      50
      command:  summarize y
 perturb vars:

----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
Perturbation results
															 Replications (50)
------------------------------------------------------------------------------
			 |   Observed        Mean          SD                  %         %
		   y |      Coef.   Perturbed   Perturbed    P-val     in CI     agree
-------------+----------------------------------------------------------------
           x |  -3.712523   -4.739402    5.343435    0.323    100.00    100.00
       _cons |  -3.374606    1.441799    5.343435    0.004     30.01     84.22
------------------------------------------------------------------------------


w tabeli:
1. Observed coef. 2. Mean 3. Diff. 4. P-val. 5. % in CI, 6. % agree
2. w e(perturb_stats): mean, median, std. dev., min, max.

*Note: Diff = observed coef. - Mean(perturbed coefficients) (niepotrzebne raczej)
Note: P-value for Ho: Diff = 0, Ha: Diff !=0
Note: "% in CI" refers to CI for observed coef.
Note: "% agree" refers to agreement to at least first significant digit

Note: "% in CI" denotes the percent of pert. coef. within the original CI
Note: "% agree" denotes the percent of pert. coef. that agree with the original
coef. to at least first significant digit
Note: P-value for Ho: Original Coef. = Mean, Ha: Original Coef. != Mean