ml_5.ado

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

*! version 3.1.5  24jan2005
program define ml_5
	version 3.1, missing
	parse "`*'", parse(" ,")
	local cmd `1'
	mac shift
	local i=length("`cmd'")
	if ("`cmd'"==substr("begin",1,`i')) { mx_begin `*' }
	else if ("`cmd'"==substr("function",1,`i')) { mx_func `*' }
	else if ("`cmd'"==substr("method",1,`i')) { mx_meth `*' }
	else if ("`cmd'"==substr("model",1,`i')) { mx_model `*' }
	else if ("`cmd'"==substr("depnames",1,`i')) { mx_depn `*' }
	else if ("`cmd'"==substr("sample",1,`i')) { mx_samp `*' }
	else if ("`cmd'"==substr("search",1,`i')) { mx_srch `*' }
	else if ("`cmd'"==substr("maximize",1,`i')) { mx_mx1 `*' }
	else if ("`cmd'"==substr("post",1,`i')) { mx_post `*' }
	else if ("`cmd'"==substr("mlout",1,`i')) { mx_out `*' }
	else if ("`cmd'"==substr("report",1,`i')) { mx_rept `*' }
	else if ("`cmd'"==substr("query",1,`i')) { mx_q `*' }
	else { error 198 }
end


*  mx_begin:  version 1.0.1  06/20/93
program define mx_begin
	version 3.1, missing 
	global S_mldepn
	global S_mlfunc
	global S_mlmeth
	global S_mlwgt
	global S_mlmb
	global S_mlsf
	global S_mlmv
	global S_mlneq
	global S_mleqn
	global S_mlsag
	global S_mlag
	global S_ll0
	global S_mlmdf
end


*  mx_d0:  version 2.0.1  25mar1996
program define mx_d0
	version 4.0, missing
	local b "`1'"       /* starting values on entrance */
	local f "`2'"       /* scalar name to contain function value */
	local g "`3'"       /* gradient vector */
	local D "`4'"       /* matrix name to contain -2nd deriv matrix */

	mac shift 4
	local options "FIRSTIT LASTIT FAST(string) *"
	parse "`*'"

	if "`lastit'" != "" { 
		capture matrix drop S_mld0S
		exit 
	}

	$S_mlfunc `b' `f'		/* f to contain f(X)	*/
	if ("`fast'"=="0") { exit }


	tempname h bb fm0 fp0 fph fmh fpp fmm

	local n = colsof(matrix(`b'))

	if "`firstit'" != "" {
		mat S_mld0S = J(1,`n',1)
	}

	matrix `h' = J(1,`n',0)		/* perturberances 	*/
	matrix `D' = J(`n',`n',0)	/* -2nd deriv matrix 	*/
	matrix `g' = `b'		/* gradients, just want names */
	matrix `fph' = J(1,`n',0)	/* will contain f(X+h_i)	*/
	matrix `fmh' = J(1,`n',0)	/* will contain f(X-h_i)	*/

	local epsf 1e-3			/* see notes below	*/

	local i 1
	while `i' <= `n' {
		mat `h'[1,`i'] = (abs(`b'[1,`i'])+`epsf')*`epsf'

		mat `bb' = `b'
		mat `bb'[1,`i'] = `b'[1,`i'] - float(S_mld0S[1,`i']*`h'[1,`i'])
		$S_mlfunc `bb' `fm0'		/* fm0 = f(X-hi)	*/

					/* adjust S, combine with h */
		adjustS0 `fm0' `f' `h' `b' `bb' `i' `firstit'


		mat `fmh'[1,`i'] = `fm0'


		mat `bb' = `b'
		mat `bb'[1,`i'] = `b'[1,`i'] + `h'[1,`i']
		$S_mlfunc `bb' `fp0'		/* fp0 = f(X+hi)	*/

		mat `fph'[1,`i'] = `fp0'	/* fph[k] = f(X+hk)	*/

						/* Gradient		*/
		mat `g'[1,`i'] = (`fp0' - `fm0')/(2*`h'[1,`i'])

						/* Dii calculation	*/
		mat `D'[`i',`i'] =-(`fp0'-2*scalar(`f')+`fm0')/((`h'[1,`i'])^2)

		local j 1
		while `j' < `i' { 
			mat `bb' = `b'
			mat `bb'[1,`i'] = `bb'[1,`i'] + `h'[1,`i']
			mat `bb'[1,`j'] = `bb'[1,`j'] + `h'[1,`j']
			$S_mlfunc `bb' `fpp'

			mat `bb' = `b'
			mat `bb'[1,`i'] = `bb'[1,`i'] - `h'[1,`i']
			mat `bb'[1,`j'] = `bb'[1,`j'] - `h'[1,`j']
			$S_mlfunc `bb' `fmm'

			mat `D'[`i',`j'] = - /*
			*/ (`fpp'+`fmm'+2*scalar(`f') /* 
			*/ -`fph'[1,`i']-`fmh'[1,`i'] /*
			*/ -`fph'[1,`j']-`fmh'[1,`j']) /*
			*/ / (2*`h'[1,`i']*`h'[1,`j'])

			mat `D'[`j',`i'] = `D'[`i',`j']
			local j=`j'+1
		}
		local i=`i'+1
	}
end

program define adjustS0 /* fm0 f h b bb i */
	version 4.0, missing
	local fm0	"`1'"	/* scalar, f(X-h[i])		*/
	local f		"`2'"	/* scalar, f(X)			*/
	local hvec	"`3'"	/* vector, value of h[]		*/
	local X		"`4'"	/* vector, base value of X[]	*/
	local XX	"`5'"	/* vector, X-h[i]		*/
	local i		"`6'"	/* i				*/

	local ep0 1e-8
	local ep1 1e-7
	local epmin 1e-10

	if "`7'"!="" {
		local ep0 1e-4
		local ep1 10e-4
	}

	local goal0 = (abs(scalar(`f'))+`ep0')*`ep0'
	local goal1 = (abs(scalar(`f'))+`ep1')*`ep1'
	local mgoal = (`goal0'+`goal1')/2
	local mingoal = (abs(scalar(`f'))+`epmin')*`epmin'

	local df = abs(scalar(`f'-`fm0'))

	tempname S newS h
	scalar `S' = S_mld0S[1,`i']
	scalar `h' = `hvec'[1,`i']
	local lastS .
	local iter 1

	while (`df'<`goal0' | `df'>`goal1') & `iter' <= 40 {
		if `df'<`mingoal' {
			scalar `newS' = `S'*2
		}
		else {
			if `lastS'>=. {
				scalar `newS'=`S'*cond(`df'<`goal0',2,.5)
			}
			else {
				/* f(0) = 0, f(S)=df, f(lastS)=lastdf */
				local a = (`lastdf'/`lastS'-`df'/`S') / /*
						*/ (`lastS'-`S')
				local b = `df'/`S' - `a'*`S'
				local r1 = ( /*
				*/ -`b' + sqrt(`b'*`b'+4*`a'*`mgoal') /* 
				*/ ) / (2*`a')
				scalar `newS'=`S'*cond(`df'<`goal0',2,.5)
				if `df'>`goal1' & `r1'>0 & `r1'<`S' {
					scalar `newS' = `r1'
				}
				else if `df'<`goal0' & `r1'>`S' & `r1'< . {
					scalar `newS' = `r1'
				}
			}
		}
		local lastS = `S'
		local lastdf= `df'
		scalar `S' = `newS'

		mat `XX' = `X'
		mat `XX'[1,`i'] = `X'[1,`i'] - float(`h'*`S')
		$S_mlfunc `XX' `fm0'		/* fm0 = f(X-hi)	*/

		local df = abs(scalar(`f'-`fm0'))
		local iter = `iter' + 1
	}
	if `df'<`goal0' | `df'>`goal1' { /* did not meet goal */
		di in red "$S_mlfunc does not compute a continuous " /*
		*/ "nonconstant function" _n /*
		*/ "could not calculate numerical derivatives"
		exit 430
	}
	matrix `hvec'[1,`i'] = float(scalar(`h')*`S')
	matrix S_mld0S[1,`i'] = `S'
end
			

/*
Derivatives and 2nd derivatives:

Notation:
	f(X)  means f(x1, x2, ...)
	f(X+h2) means f(x1, x2+h2, x3, ...)
	f(x+h_2+h_3) means f(x1, x2+h2, x3+h3,...)  etc.

Centered first derivatives are:

	g_i = [f(X+hi)-f(X-hi)]/2hi

Centered d^2 f/dx_i^2 are

	Dii = [ f(X+hi) - 2*f(X) + f(X-hi) ] / (hi)^2

Cross partials are:

	Dij = [ f(++) + f(--) + 2*f(00) - f(+0) - f(-0) - f(0+) - f(0-) ]
			/ (2*hi*hj)

We define hi = (|xi|+eps)*eps where eps is sqrt(machine precision) or
maybe cube root.

Program logic:
	calculate f(X)
	for i { 
		calculate hi
		calculate f(X-hi), optimizing hi, and save 
		calculate f(X+hi) and save
		obtain g_i
		calculate Dii
		for j<i { 
			calculate Dij and save in Dji too
		}
	}




For an even more expensive calculation of the 2nd, using formula:

	Dij = [f(++) - f(+-) - f(-+) + f(--)]/(4*hi*hj)

use the code:

			mat `bb' = `b'
			mat `bb'[1,`i'] = `bb'[1,`i'] + `h'[1,`i']
			mat `bb'[1,`j'] = `bb'[1,`j'] + `h'[1,`j']
			$S_mlfunc `bb' `fpp'

			mat `bb' = `b'
			mat `bb'[1,`i'] = `bb'[1,`i'] + `h'[1,`i']
			mat `bb'[1,`j'] = `bb'[1,`j'] - `h'[1,`j']
			$S_mlfunc `bb' `fpm'

			mat `bb' = `b'
			mat `bb'[1,`i'] = `bb'[1,`i'] - `h'[1,`i']
			mat `bb'[1,`j'] = `bb'[1,`j'] + `h'[1,`j']
			$S_mlfunc `bb' `fmp'

			mat `bb' = `b'
			mat `bb'[1,`i'] = `bb'[1,`i'] - `h'[1,`i']
			mat `bb'[1,`j'] = `bb'[1,`j'] - `h'[1,`j']
			$S_mlfunc `bb' `fmm'

			mat `D'[`i',`j'] = - /*
				*/ (`fpp'-`fpm'-`fmp'+`fmm')/ /*
				*/ (4*`h'[1,`i']*`h'[1,`j'])
*/


*  mx_d1:  version 1.2.0  01may1995
program define mx_d1
	version 4.0, missing
	local b "`1'"       /* starting values on entrance */
	local d "`3'"       /* derivative matrix */
	local f "`2'"       /* scalar name to contain function value */
	local v "`4'"       /* matrix name to contain variance matrix */

	mac shift 4
	local options "LASTIT FIRSTIT FAST(string)"
	parse "`*'"

	if ("`lastit'"~="") { exit }

	$S_mlfunc `b' `f' `d' , fast(`fast')
	if (`fast'==0) { exit }

	tempname bb dd1 ff h dd2
	local epsf 1e-4			/* was 1e-3, ought to be settable */
	cap mat drop `v'
	local i 1
	while (`i' <= colsof(matrix(`b'))) {
		scalar `h' = float((abs(`b'[1,`i'])+`epsf')*`epsf')
		mat `bb' = `b'
		mat `bb'[1,`i'] = `b'[1,`i'] + `h'
		$S_mlfunc `bb' `ff' `dd1'

		mat `bb' = `b'
		mat `bb'[1,`i'] = `b'[1,`i'] - `h'
		$S_mlfunc `bb' `ff' `dd2'

		scalar `h' = 1/(2*`h')
		mat    `dd2' = `dd2' - `dd1'
		mat    `dd2' = `h' * `dd2'
		mat    `v' = `v' \ `dd2'

		local i = `i' + 1
	}
	/* symmetrize */
	mat        `dd2' = `v' '
	mat        `v' = `v' + `dd2'
	mat        `v' = `v' * 0.5

	if "$S_mldbug" != "1" { exit }

/* debug code for user ml functions */
	mat colnames `d' = `pnames'
	mat coleq    `d' = `peq'
	mat colnames `v' = `pnames'
	mat coleq    `v' = `peq'
	mat rownames `v' = `pnames'
	mat roweq    `v' = `peq'
	mat list `d'
	mat list `v'
end


*  mx_depn:  version 1.0.1  06/20/93
program define mx_depn
	version 3.1, missing
	local varlist "req ex"
	parse "`*'"
	global S_mldepn "`varlist'"
end



*  mx_func:  version 1.0.1  07/20/93
program define mx_func
	version 3.1, missing
	if ("`1'"=="" | "`2'"!="") {
		error 198
	}
	global S_mlfunc `1'
end



*  mx_lf:  version 2.0.2  25mar1996
program define mx_lf
	version 4.0, missing
	local b "`1'"       /* starting values on entrance */
	local d "`3'"       /* derivative vector */
	local f "`2'"       /* matrix name to contain function value */
	local h "`4'"       /* matrix name to contain hessian matrix */

	mac shift 4
	local options "FAST(string) LASTIT FIRSTIT"
	parse "`*'"

	if "`firstit'" != "" {
		setup `b'
	}
	else if "`lastit'" != "" {
		alldone
		exit
	}

	tempname wrk
	tempvar ff sumff

	local i 1
	while `i'<=$S_mlfeq {
		tempname x`i'
		local se : word `i' of $S_mlfsd
		local ee : word `i' of $S_mlfed
		mat `wrk' = `b'[1,`se'..`ee']
		mat score double `x`i'' = `wrk' if $S_mlwgt
		local arglist "`arglist' `x`i''"
		local i = `i' + 1
	}

	qui gen double `ff' = .
	$S_mlfunc `ff' `arglist'
	qui count if `ff'>=. & $S_mlwgt!=0
	if _result(1)!=0 {
		scalar `f' = -1e30
		exit /* impossible value exit */
	}
	gen double `sumff' = sum(`ff'*$S_mlwgt)
	scalar `f' = `sumff'[_N]
	if `fast'==0 { exit }

	drop `sumff'
				/* we now continue to make derivative
				   calculations
				*/
	local epsf 1e-4
	tempname nfac
	tempvar one x0 grad xj0 fphh fmhh Dii Dij

				/* set initial scale to x`i' */
	if "`firstit'" != "" { 
		local i 1 
		while `i'<=$S_mlfeq { 
			local S : word `i' of $S_mlfS
			scalar `S' = 1
			local i=`i'+1
		}
	}

	mat `d' = J(1,$S_mlfn,0)
	mat `h' = J($S_mlfn,$S_mlfn,0)
	qui gen byte `one' = 1 if $S_mlwgt

	quietly {
		local i 1
		while `i'<=$S_mlfeq {
/* using wrk 	*/
			local se : word `i' of $S_mlfsd
			local ee : word `i' of $S_mlfed
			mat `wrk' = `b'[1,`se'..`ee']
			local vl`i' : colnames(`wrk')
			local i_cons = colnumb(`wrk',"_cons") 
			if `i_cons'< .  {
				parse "`vl`i''", parse(" ")
				local `i_cons' "`one'"
				local vl`i' "`*'"
			}
/* wrk now free */

			tempname h`i'
			tempvar fph`i' fmh`i'

			local S : word `i' of $S_mlfS
			qui summ `x`i'' [aw=$S_mlwgt]
			scalar `h`i'' = float((abs(_result(3))+`epsf')*`epsf')

			rename `x`i'' `x0'
			gen double `x`i'' = `x0' - float(`h`i''*`S')
			gen double `fmh`i'' = .