ml_opt.ado

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

*! version 7.2.20  09feb2005
program define ml_opt
	version 6 
	syntax [, DIFficult EIGEN(real 1e-7) noWARNing]
		/* eigen() should not have to be reset; it only applies
		   to -difficult- method.
		*/

	local sayback 6  /* print out "(backed up)" if it backed up `sayback' or
	                    more times */

	if "`difficu'"=="" {
		local inverse "Inverse"
		local ncstep  "NCStep"
	}
	else {
		local inverse "DiffInv"
		local ncstep  "DiffStep"
	}

	if $ML_trace==1       { local dilike "DiLike" }
	else if $ML_trace==0  { local dilike "*" }

	if $ML_trace>1 | $ML_dider  { local dilike "DiLikeTr" }
	else local tr "*"
		/* `tr'=="" => display BIG iteration log, set off by dashes */

	if "$ML_nrsho" != "" & $ML_trace != 0 {
		local dinrval "DiNRvalue"
	}
	else	local dinrval "*"

	if $ML_trace==2 | $ML_trace==4 { local dicoef "DiCoef" }
		/* display coefficient vector */

	if $ML_dider | $ML_trace > 2 { local dider "DiDer" }
	else local dider "*"
		/* display gradient length, gradient vector, and/or hessian */

	if $ML_dider < 2 { local saveH "*" }

	`tr'   di in smcl in gr "{hline 78}" _n "Iteration 0:"
	`dicoef'

	if $ML_C {	/* initial points for constraint model */ 
		mat $ML_b = $ML_b*$ML_CT
		mat $ML_b = $ML_b*$ML_CT' + $ML_Ca
	}
		
	$ML_eval 2
	if scalar($ML_f)==. {
		di in red "initial values not feasible"
		exit 1400
	}
	ml_log
	
	if $ML_C { 		/* unconstraint => constraint */
		mat $ML_b = ($ML_b-$ML_Ca)*$ML_CT
		mat $ML_V = $ML_CT'*$ML_V*$ML_CT
		mat $ML_g = $ML_g*$ML_CT
	}


	tempname b0 f0 H gPg bPb
	local kback 0
	local nc_ct 0
	local back_ct 0
	local ic 1
	local conv 0
	while !`conv' & `ic'<=$ML_iter {
		scalar `f0' = scalar($ML_f)
		mat `b0' = $ML_b
		`saveH' mat `H' = $ML_V
		`inverse'
		local notcc = r(notcc)
		if `notcc' {
			`dilike' "(not concave)"
			`dinrval' `nrval'
			`dider' `H'
			`ncstep' `eigen'
			local kback 0
			local nc_ct = `nc_ct' + 1
		}
		else {
			if `kback' >= `sayback' {
				`dilike' "(backed up)"
				`dinrval' `nrval'
					/* when it stepped to this point, it
					   did `sayback' or more step halvings
					*/
				local back_ct = `back_ct' + 1
			}
			else {
				`dilike'
				`dinrval' `nrval'
				local back_ct 0
			}

			`dider' `H'
			Step
			local kback = r(k_back)
			local nc_ct = 0
		}

		`tr' di in smcl in gr "{hline 78}" _n "Iteration `ic':"
		`dicoef'

		if `ic' == 1 & $ML_iton1 {
			continue, break
		}
		CnsEval 2

		if scalar($ML_f) == . {
			di in red "discontinuous region encountered"
			di in red "cannot compute an improvement"
			exit 430
		}

		if `ic' == 1 {
			chk4scores
		}

		ml_log

		local conv = mreldif(matrix($ML_b),`b0')<=$ML_tol | /*
		*/ reldif(scalar($ML_f),`f0')<=$ML_ltol

		if `conv' {
			if "$ML_gtol" != "" {
				ChkGrad
				local conv = ( r(grad_ok) /*
				*/ | `back_ct' > 30 | `nc_ct' > 30)
			}
		}
		if (`conv' & "$ML_nrtol" != "") {
			ChkNRtol conv nrval :
			if `conv' {
				local kback 0
			}
		}
		else {
			local nrval
		}

		local ic = `ic' + 1
	}

	`saveH' mat `H' = $ML_V
	if "$ML_noinv" == "" {
		capture mat $ML_V = syminv($ML_V)
		if _rc {
			di in red "Hessian has become unstable or asymmetric"
			exit _rc
		}
	}
	if diag0cnt(matrix($ML_V)) {
		`dilike' "(not concave)"
		`dinrval' `nrval'
	}
	else if `kback' >= `sayback' {
		`dilike'  "(backed up)"
		`dinrval' `nrval'
	}
	else {
		`dilike'
		`dinrval' `nrval'
	}
	`dider' `H'

	`tr' di in smcl in gr "{hline 78}"

	if `conv' {
		global ML_rc 0
		global ML_conv 1
	}
	else if "`warning'" == "" {
		if $ML_rc { di in red "convergence not achieved" }
		else	  { di in blu "convergence not achieved" }
	}
	
	if $ML_C { 		/* constraint => unconstraint */
		mat $ML_b = $ML_b*$ML_CT' + $ML_Ca
		mat $ML_V = $ML_CT*$ML_V*$ML_CT'
		mat $ML_g = $ML_g*$ML_CT'
	}

end

program chk4scores
	version 8.0
	if "$ML_sclst" != "" & ///
	("$ML_vce" == "robust" | "$ML_vce2" == "OPG") {
		local opt = cond("$ML_vce"=="robust","robust","vce(opg)")
		local ok 0
		foreach var of global ML_sclst {
			capture assert missing(`var')
			if (_rc) local ++ok
		}
		if !`ok' {
			di as err ///
	"$ML_user failed to compute scores required by the `opt' option"
			exit 504
		}
	}
end

program define DiLike /* "message" */
	local ic = $ML_ic - 1
	di in smcl in gr "Iteration `ic':" _col(16) "$ML_crtyp = " /*
	*/ in ye %10.0g scalar($ML_f) in blu "  `*'"
end

program DiNRvalue
	args val
	if "`val'" == "" {
		exit
	}
	local nrname "g inv(H) g'"
	local len1 = length("$ML_crtype")
	local len2 = length("`nrname'")
	if $ML_trace == 1 {
		local col = 16 + max(`len1',`len2') - 1
		di in smcl as txt _col(16) "`nrname'" _col(`col') " = " ///
			in ye %10.0g `val'
	}
	else {
		local col = 66 - `len2'
		di in smcl as txt _col(`col') "`nrname' = " ///
			in ye %10.0g `val'
	}
end

program define DiLikeTr /* "message" */
	*if $ML_trace == 2 { local nl "_n" }
	local col = 66-length("$ML_crtyp")
	di `nl' in gr _col(`col') "$ML_crtyp = " in ye %10.0g scalar($ML_f)

	if "`*'"!="" {
		local col = 79 - length("`*'")
		di in blu _col(`col') "`*'"
	}
end

program define DiCoef
	di in gr "Coefficient vector:"
	mat list $ML_b, noheader noblank format(%9.0g)
	di /* blank line */
end

program define DiDer
	args H
	tempname c
	mat `c' = $ML_g*$ML_g'

	if $ML_dider == 0 | ($ML_dider == 2 & $ML_trace == 3) {
		di in gr "Gradient vector length =" in ye %9.0g sqrt(`c'[1,1])
		if $ML_dider == 0 { exit }
		di /* blank line */
	}

	if $ML_dider == 1 {
		_cpmatnm $ML_b, vec($ML_g)
	}
	else if $ML_dider == 2 {
		_cpmatnm $ML_b, square(`H')
	}
	else	_cpmatnm $ML_b, vec($ML_g) square(`H')

	if ($ML_dider == 1 | $ML_dider == 3) /*
	*/ & (substr("$ML_meth",3,.)!="debug" | $ML_trace == 4) {
		di in gr "Gradient vector (length =" in ye %9.0g /*
		*/ sqrt(`c'[1,1]) in gr "):"
		mat list $ML_g, noheader noblank format(%9.0g)
		local newline "_n"
	}
	if $ML_dider > 1 & substr("$ML_meth",3,.)!="debug" {
		di `newline' in gr "Negative Hessian:"
		mat list `H', noheader noblank format(%9.0g)
	}
end

program define Inverse, rclass
	tempname H
	mat rename $ML_V `H', replace
	if "$ML_noinv" == "" {
		capture mat $ML_V = syminv(`H')
		if _rc {
			di in red "Hessian has become unstable or asymmetric"
			exit _rc
		}
	}
	else	mat $ML_V = `H'
	IsOk `H' $ML_V
	if r(okay) {
		return scalar notcc = 0	/* meaning not concave is FALSE */
		exit
	}
	return scalar notcc = 1		/* meaning not concave is TRUE	*/
	if "$ML_noinv" != "" {
		capture mat `H' = syminv($ML_V)
		if _rc {
			di in red "Hessian has become unstable or asymmetric"
			exit _rc
		}
	}

	local origt = trace($ML_V)
	local amt 1.1			/* amt must be > 1	*/
	while !r(okay) {
		local dim = matrix(colsof(`H'))
		local i 1
		while `i' <= `dim' {
			mat `H'[`i',`i'] = `H'[`i',`i']+`amt'*abs(`H'[`i',`i'])
			local i = `i' + 1
		}
		local amt = 2*`amt'
		capture mat $ML_V = syminv(`H')
		if _rc {
			di in red "Hessian has become unstable or "	/*
				*/ "asymmetric (NC)"
			exit _rc
		}
		IsOk `H' $ML_V
	}
	local adj = `origt'/trace($ML_V)
	if `adj' > 0 {
		mat $ML_V = `adj'*$ML_V
	}
end

program define DiffInv, rclass
	tempname Hinv
	if "$ML_noinv" == "" {
		capture mat `Hinv' = syminv($ML_V)
		if _rc {
			di in red "Hessian has become unstable or "	/*
				*/ "asymmetric (D)"
			exit _rc
		}
	}
	else	mat `Hinv' = $ML_V
	IsOk $ML_V `Hinv'
	if r(okay) {
		mat rename `Hinv' $ML_V, replace
		return scalar notcc = 0	/* meaning not concave is FALSE */
		exit
	}
	return scalar notcc = 1		/* meaning not concave is TRUE	*/
	if "$ML_noinv" != "" {
		capture mat `H' = syminv($ML_V)
		if _rc {
			di in red "Hessian has become unstable or asymmetric"
			exit _rc
		}
	}
end

program define IsOk /* H invH */, rclass
	args H V
	if diag0cnt(matrix(`V'))==0 {
		return scalar okay = 1
		exit
	}
	if diag0cnt(matrix(`H'))==0 {
		return scalar okay = 0
		exit
	}
	local dim = matrix(colsof(`H'))
	local i 1
	while `i' <= `dim' {
		if `V'[`i',`i']==0 & `H'[`i',`i']!=0 {
			return scalar okay = 0
			exit
		}
		local i = `i' + 1
	}
	return scalar okay = 1
end

program define Step, rclass
	if $ML_trace < 3 { local show "*" }

	tempname f0 b0 d
	scalar `f0' = scalar($ML_f)
	mat `b0' = $ML_b
	mat `d' = $ML_g*$ML_V
	mat $ML_b = `d' + $ML_b

	`show' DiStep `d'

	CnsEval 0
	`show' DiLikeTr

	if scalar($ML_f)==. | scalar($ML_f)<`f0' {
		`show' di in blu _col(61) "(initial step bad)"
		Backward `f0' `b0' 60
		return scalar k_back = r(k_back)
		exit
	}

	`show' di in blu _col(60) "(initial step good)"

	Forward `f0' `b0' `d' 0.125 2

	return scalar k_back = 0
end

program define DiStep
	args step
	tempname c
	mat `c' = `step'*`step''

	if $ML_trace == 3 {
		if !$ML_dider {
			di in gr "Step length            =" /*
			*/ in ye %9.0g sqrt(`c'[1,1]) _n /*
			*/ in gr "Stepping b + step -> new b"
		}
		else {
			di _n in gr "Step length =" /*
			*/ in ye %9.0g sqrt(`c'[1,1]) _n /*
			*/ in gr "Stepping b + step -> new b"
		}
		exit
	}

	_cpmatnm $ML_b, vec(`step')
	di _n in gr "Step (length =" in ye %9.0g sqrt(`c'[1,1]) in gr "):"
	mat list `step', noheader noblank format(%9.0g)

	di _n in gr "b + step -> new b:"
	mat list $ML_b, noheader noblank format(%9.0g)
	di /* blank line */
end

program define NCStep, rclass
	if $ML_trace < 3 { local show "*" }

	tempname f0 b0 d
	scalar `f0' = scalar($ML_f)
	mat `b0' = $ML_b
	mat `d' = $ML_g*$ML_V    /* Newton-Raphson step */

	Stepsize `f0' $ML_g `d'  /* `d' = scaled Newton-Raphson step. */

	mat $ML_b = `d' + $ML_b

	`show' DiStep `d'

	CnsEval 0