permute.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

	char _dta[seed] "`seed'"
	char _dta[permvar] "`permvar'"
	char _dta[command] "`command'"
	char _dta[pm_version] 2
	quietly drop in 1

	if `"`filetmp'"' == "" {
		quietly save `"`saving'"', replace
	}

	ClearE
	Results, `ropts'
end

program Results
	syntax [anything(name=namelist)]	///
		[using/] [,			///
			eps(real 1e-7)		/// -GetResults- options
			left			///
			right			///
			TItle(passthru)		///
			Level(cilevel)		/// -DisplayResults- options
			noHeader		///
			noLegend		///
			Verbose			///
			notable			/// not documented
		]

	if `"`using'"' != "" {
		preserve
		qui use `"`using'"', clear
	}
	else if `"`namelist'"' != "" {
		local namelist : list uniq namelist
		preserve
	}
	local 0 `namelist'
	syntax [varlist(numeric)]
	if "`namelist'" != "" {
		keep `namelist'
		local 0
		syntax [varlist]
	}

	GetResults `varlist',	///
		eps(`eps')	///
		`left' `right'	///
		level(`level')	///
		`title'
	DisplayResults, `header' `table' `legend' `verbose'
end

program GetResults, rclass
	syntax varlist [,		///
		Level(cilevel)		///
		eps(real 1e-7)		///
		left			///
		right			///
		TItle(string asis)	///
	]

	// get data characteristics
	// data version
	local version : char _dta[pm_version]
	capture confirm integer number `version'
	if c(rc) | "`version'" == "" {
		local version 1
	}
	else if `version' <= 0 {
		local version 1
	}
	// original number of observations
	local obs : char _dta[N]
	if "`obs'" != "" {
		capture confirm integer number `obs'
		if c(rc) {
			local obs
		}
		else if `obs' <= 0 {
			local obs
		}
	}
	// number of strata
	local nstrata : char _dta[N_strata]
	if "`nstrata'" != "" {
		capture confirm integer number `nstrata'
		if c(rc) {
			local nstrata
		}
	}
	// strata variable
	if "`nstrata'" != "" {
		local strata : char _dta[strata]
		if "`strata'" != "" {
			capture confirm names `strata'
			if c(rc) {
				local strata
			}
		}
	}
	// permutation variable
	local permvar : char _dta[permvar]
	capture confirm name `permvar'
	if c(rc) | `:word count `permvar'' != 1 {
		local permvar
	}
	if `"`permvar'"' == "" {
		di as error ///
"permutation variable name not present as data characteristic"
		exit 9
	}

	// requested event
	GetEvent, `left' `right' eps(`eps')
	local event `s(event)'
	local rel `s(rel)'
	local abs `s(abs)'

	tempvar diff
	gen `diff' = 0
	local K : word count `varlist'
	tempname b c reps p se ci
	matrix `b' = J(1,`K',0)
	matrix `c' = J(1,`K',0)
	matrix `reps' = J(1,`K',0)
	matrix `p' = J(1,`K',0)
	matrix `se' = J(1,`K',0)
	matrix `ci' = J(1,`K',0) \ J(1,`K',0)

	local seed : char _dta[seed]
	local k_eexp 0
	forvalues j = 1/`K' {
		local name : word `j' of `varlist'
		local value : char `name'[permute]
		capture matrix `b'[1,`j'] = `value'
		if c(rc) | missing(`value') {
			di as err ///
`"estimates of observed statistic for `name' not found"'
			exit 111
		}
		quietly replace ///
		`diff' = (`abs'(`name') `rel' `abs'(`value') - `eps')
		sum `diff' if `name'<., meanonly
		if r(N) < c(N) {
			local missing missing
		}
		mat `c'[1,`j'] = r(sum)
		mat `reps'[1,`j'] = r(N)
		quietly cii `=`reps'[1,`j']' `=`c'[1,`j']', level(`level')
		mat `p'[1,`j'] = r(mean)
		mat `se'[1,`j'] = r(se)
		mat `ci'[1,`j'] = r(lb)
		mat `ci'[2,`j'] = r(ub)
		local coleq `"`coleq' `"`:char `name'[coleq]'"'"'
		local colname `colname' `:char `name'[colname]'
		if `version' >= 2 {
			local exp`j' : char `name'[expression]
		}
		if `"`:char `name'[is_eexp]'"' == "1" {
			local ++k_eexp	
		}
	}
	local coleq : list clean coleq

	if `version' >= 2 {
		// command executed for each permutation
		local command : char _dta[command]
		local k_exp = `K' - `k_eexp'
	}
	else {
		local k_eexp 0
		local k_exp 0
	}

	// put stripes on matrices
	if `"`coleq'"' == "" {
		matrix colnames `b' = `varlist'
	}
	else {
		matrix colnames `b' = `colname'
		if `"`coleq'"' != "" {
			matrix coleq `b' = `coleq'
		}
	}
	matrix rowname `b' = y1
	_copy_mat_stripes `c' `reps' `p' `se' `ci' : `b', novar
	matrix rowname `ci' = ll ul
	matrix roweq `ci' = _ _

	// Save results
	return clear
	return scalar version = `version'
	if "`obs'" != "" {
		return scalar N = `obs'
	}
	return scalar level = `level'
	return scalar k_eexp = `k_eexp'
	return scalar k_exp = `k_exp'
	return matrix reps `reps'
	return matrix c `c'
	return matrix b `b'
	return matrix p `p'
	return matrix se `se'
	return matrix ci `ci'
	return local seed `seed'
	return local missing `missing'
	return local permvar `permvar'
	if "`nstrata'" != "" {
		return scalar N_strata = `nstrata'
		if "`strata'" != "" {
			return local strata `strata'
		}
	}
	return local event `event'
	return local left `left'
	return local right `right'
	forval i = 1/`K' {
		return local exp`i' `"`exp`i''"'
	}
	if `"`title'"' != "" {
		return local title `"`title'"'
	}
	else	return local title "Monte Carlo permutation results"
	return local command `"`command'"'
	return local cmd permute
end

program DisplayResults, rclass
	syntax [,			///
		noHeader		///
		noLegend		///
		Verbose			///
		notable			///
	]

	if "`header'" == "" {
		_coef_table_header, rclass
		if r(version) >= 2 & "`legend'" == "" {
			_prefix_legend permute, rclass `verbose'
			di as txt %`s(col1)'s "permute var" ":  `r(permvar)'"
		}
	}

	// NOTE: _coef_table_header needs the results in r() to work properly,
	// thus the following line happens here instead of at the very top.
	return add

	if ("`table'" != "") {
		exit
	}
	else if "`header'" == "" {
		di
	}

	tempname Tab results
	.`Tab' = ._tab.new, col(8) lmargin(0) ignore(.b)
	ret list
	// column           1      2     3     4     5     6     7     8
	.`Tab'.width	   13    |12     8     8     8     8    10    10
	.`Tab'.titlefmt %-12s      .     .     .     .     .  %20s     .
	.`Tab'.pad	    .      2     0     0     0     0     0     1
	.`Tab'.numfmt       .  %9.0g     .     . %7.4f %7.4f     .     .

	local cil `=string(`return(level)')'
	local cil `=length("`cil'")'
	if `cil' == 2 {
		local cititle "Conf. Interval"
	}
	else {
		local cititle "Conf. Int."
	}
                                                                                
	// begin display
	.`Tab'.sep, top
	.`Tab'.titles "T" "T(obs)" "c" "n" "p=c/n" "SE(p)" ///
		"[`return(level)'% `cititle']" ""

	tempname b c reps p se ci
	matrix `b' = return(b)
	matrix `c' = return(c)
	matrix `reps' = return(reps)
	matrix `p' = return(p)
	matrix `se' = return(se)
	matrix `ci' = return(ci)
	local K = colsof(`b')
	local colname : colname `b'
	local coleq   : coleq `b', quote
	local coleq   : list clean coleq
	if `"`:list uniq coleq'"' == "_" {
		local coleq
	}
	gettoken start : colname
	forvalues j = 1/`K' {
		local name : word `j' of `colname'
		local curreq : word `j' of `coleq'
		if "`curreq'" != "`eq'" {
			.`Tab'.sep
			di as res %-12s abbrev("`curreq'",12) as txt " {c |}"
			local eq `curreq'
		}
		else if "`name'" == "`start'" {
			.`Tab'.sep
		}
		local name = abbrev(`"`name'"',10)
		.`Tab'.row "`name'"	///
			`b'[1,`j']	///
			`c'[1,`j']	///
			`reps'[1,`j']	///
			`p'[1,`j']	///
			`se'[1,`j']	///
			`ci'[1,`j']	///
			`ci'[2,`j']	///
			// blank
	}
	.`Tab'.sep, bottom
	TableFoot "`return(event)'" `K' `return(missing)'
end

program ClearE, eclass
	ereturn clear
end

program GetEvent, sclass
	sret clear
	syntax [, left right eps(string)]
	if "`left'"!="" & "`right'"!="" {
		di as err "only one of left or right can be specified"
		exit 198
	}
	if "`left'"!="" {
		sreturn local event "T <= T(obs)"
		sreturn local rel "<="
		sreturn local eps = -`eps'
	}
	else if "`right'"!="" {
		sreturn local event "T >= T(obs)"
		sreturn local rel ">="
	}
	else {
		sreturn local event "|T| >= |T(obs)|"
		sreturn local rel ">="
		sreturn local abs "abs"
	}
end

program PermVars // "byvars" k var
	args strata k x
	tempvar r y
	quietly {
		if "`strata'"!="" {
			by `strata': gen double `r' = uniform()
		}
		else	gen double `r' = uniform()

		sort `strata' `r'
		local type : type `x'
		gen `type' `y' = `x'[`k']
		drop `x'
		rename `y' `x'
	}
end

program PermDiV // "byvars" k min max var
	args strata k min max x
	tempvar y
	if "`strata'"!="" {
		sort `strata'
		local bystrata "by `strata':"
	}
	quietly {
		gen byte `y' = . in 1
		`bystrata' replace `y' = ///
			uniform()<(`k'-sum(`y'[_n-1]))/(_N-_n+1)
		replace `x' = cond(`y',`max',`min')
	}
end

program TableFoot
	args event K missing
	if `K' == 1 {
		di as txt ///
"Note:  Confidence interval is with respect to p=c/n."
	}
	else {
		di as txt ///
"Note:  Confidence intervals are with respect to p=c/n."
	}
	if "`event'"!="" {
		di in smcl as txt "Note:  c = #{`event'}"
	}
	if "`missing'" == "missing" {
		di as txt ///
"Note:  Missing values observed in permutation replicates."
	}
end

exit

----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
Logit estimates                                   Number of obs   =        100
                                                  LR chi2(1)      =       0.83
                                                  Prob > chi2     =     0.3610
Log likelihood = -5.1830091                       Pseudo R2       =     0.0745
        
------------------------------------------------------------------------------
           y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |  -3.712523   4.739402    -0.78   0.433    -13.00158    5.576535
       _cons |  -3.374606   1.441799    -2.34   0.019    -6.200481   -.5487309
------------------------------------------------------------------------------
        
----+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8
Monte Carlo Permutation Statistics                Number of obs    =        74
                                                  Replications     =       100

------------------------------------------------------------------------------
T            |     T(obs)       c       n   p=c/n   SE(p) [95% Conf. Interval]
-------------+----------------------------------------------------------------
_pm1         |   24.77273       0     100  0.0000  0.0000         0   .0362167 
------------------------------------------------------------------------------
Note:  confidence interval(s) are with respect to p=c/n
Note:  c = #{|T| >= |T(obs)|}